AD-A163  572  IMPROVED  NUMERICAL  IMPLEMENTATION  OF  THE  BOUNDING  1/i 

SURFACE  PLASTICITV  MODEL  FOR  COHESIVE  SOILS(U) 

'  CALIFORNIA  UNIV  DAVIS  L  R  HERRMANN  ET  AL  DEC  85 

UNCLASSIFIED  NCEL-CR-8S  004  N62583-85-M-T176  F/G  8/13  NL 


NCEL 

Contract  Report 


CR-86.004 


December  1985 
An  Investigation  Conducted  By 
University  Of  California  Davis.  California 

Sponsored  By  Naval  Facilities 
Engineering  Command 


IMPROVED  NUMERICAL 
IMPLEMENTATION  OF  THE 
BOUNDING  SURFACE 
PLASTICITY  MODEL 
FOR  COHESIVE  SOILS 


ABSTRACT  A  substantially  more  robust  numerical  algorithm  for  the 
evaluation  of  the  bounding  surface  plasticity  model  for  cohesive  soils  is 
developed.  The  robustness  of  the  new  algorithm  assures  accurate  results 
for  reasonable  sized  solution  steps  and  qualitatively  correct  predictions, 
even  for  exceptionally  large  steps.  The  improved  predictions  are 
illustrated  for  three  oxamplns. 


t  L  *? 

\  JAN  b  i 


NAVAL  CIVIL  ENGINEERING  LABORATORY  PORT  HUENEME  CALIFORNIA  93043 


TbI*  do::’  . 
fc»  pu; 

dintilL  .  ;i 


1  31  0  2  0 


METRIC  CONVERSION  FACTORS 


SirVS  "A'S  3*  * 


2  8#  ill  s*i  I  1 1 

sl  lliti  m§  !li  iilli?  *i 


fX|  S  «  n  »  a  S  n  «  ia  4  §  r«  iu|5  —  8  8  n  UJ  |s 
h  bon-d  <|d-d«  •  o  m  -  J  d«-d«i-  K  «>  g 

I  i  s  i  3  *2  s? 


HIS 


>  5S  i  5  5  -i  -  «  -  “  ;  5 

e  Si  5  Bcjtw  E  —  i  88  k 

1  iiiii  fill  HI  mm  It 


I  §  e  e  J  *WSl  .5-  ? _ Ve  o° 


02  1 81  181  111  191  91  |»l  I  Cl  |21  |U  101  IB  18 


I  I  I  I  I  I  I  I  I 


111  L 


TjfP  'Vjr  *V|T  Tf  V  TfP  T|‘lT|,|TP1‘|‘l,|Tp’rr|T|’l,|,l‘  T]TPP|T  T|TfP|  ■  1 1 1  ■  1 1 


§  §  E I  Te^eI^  .5.  11? _ ”e"e  o° 


I  S 

5  2  5  i 

•3  •  •  C 


1 1  if  |  ml.  i 

'  Hu  mil  in 

»  2 

X  Si  Ui 

•1  •£!  2  *•  o><o“io8«.»i  $  a  3 

I  I  o 


5  5  2 

III  8 1  3 1 

llllllla j  2|3S 


IAJ  ? 

ro  to  *  S  -3 

ts  a  3  .£  u 

M  A  «  H  <*  • 


J  ifliaaaoonoo  5  c  2 
>1  X 


t  In!  .if  s|  s|  If 

■  lit?  mil  mi  ijiiiu??  if 


sell  'VrVl  32 


SECURITY  CLASSIFICATION  OF  This  RAGE  C**#*  Hot o  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1  REPORT  NUMBER  J  GOVT  ACCESSION  NO. 

CR  86.004  Ah.AUl 

3  RECIPIENT’S  CATALOG  NUMBER 

A.  TITLE  (and  Subtitle) 

Improved  Numerical  Implementation  of 
the  Bounding  Surface  Plasticity  Model 
for  Cohesive  Soils 

S  type  of  REPORTjt  PERIOO  COVEREO 

Jan  1985°-  Sep  1985 

«  PERFORMING  ORG.  REPORT  NUMBER 

’  Leonard  R.  Herrmann,  PhD 

Victor  Kaliakin 

C.K.  Shen,  PhD 

a  CONTRACT  or  grant  NUMBERYtJ 

N62583/85MT176 

'i  ;\ 

»  PERFORMING  ORGANIZATION  NAME  ano  aooress 

University  of  California 

Davis,  CA 

to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  &  WORK  UNIT  NUMBERS 

61153N 

YR023.03.01.006 

It.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Naval  Civil  Engineering  Laboratory 

Port  Hueneme,  CA  93043-5003 

>2.  REPORT  DATE 

December  1985 

n.  NUMBER  OF  PAGES 

64 

14.  MONITORING  AGENCY  name  4  AOORESS^/  different  front  Controlling  Office) 

Naval  Facilities  Engineering  Command 

200  Stovall  Street 

Alexandria,  VA  22332-2300 

IS-  SECURITY  CLASS,  (o  1  th to  roport ) 

Unclassi fied 

1S«.  OECL  ASSIFIC  ATtON  DOWNGRADING 
SCHEDULE 

IS.  DISTRIBUTION  STATEMENT  (of  thie  Report) 

Approved  for  public  release;  distribution  is  unlimited. 

17.  DISTRIBUTION  STATEMENT  (of  the  ebetract  entered  in  Block  20.  •!  different  from  Roport) 

is.  supplementary  notes 

^finite  elements  constitutive  lawf  soil  plasticity;  cohesive 
soil;  numerical  implementation,^ 

f 

20^  ABSTRACT  ( Conimua  on  roeere*  •  ido  II  necettary  end  identity  by  block  number) 

/A  substantially  more  robust  numerical  algorithm  for  the 
evaluation  of  the  bounding  surface  plasticity  model  for 
cohesive  soils  is  developed.  The  robustness  of  the  new 
algorithm  assures  accurate  results  for  reasonably  sized  solution 
steps  and  qualitatively  correct  predictions,  even  for  exception¬ 
ally  large  steps.  The  improved  predictions  are  illustrated  for 
t.hrpp  pyamnlpc  . 

oo  , 1473  eotrioH of. i movm.j obsolete  Unclassified 

s  E  CuRITY  CL  ASSlFlCATlON-OF  THIS  PAGE  Dato  Enf  er  ed) 


'«JU1  IJJJ.  ^.1, 


CONTENTS 

Page 


INTRODUCTION  .  1 

PRELIMINARIES  .  2 

MODIFICATIONS  TO  "CLAY"  .  5 

IMPLEMENTATION  OF  CLAY .  13 

EXAMPLES .  24 

CONCLUSIONS .  25 

REFERENCES .  31 

APPENDIX  I  .  34 

APPENDIX  II  .  51 

APPENDIX  III .  55 

LIST  OF  FIGURES 

Figure  1.  Schematic  illustration  of  a  bounding  surface  in 

general  stress  space  .  26 

Figure  2.  Schematic  illustration  of  a  cross-section  of  the 
bounding  surface  for  cohesive  soils  in  invariant 
stress  space  .  27 

Figure  3.  Plane  strain  footing  example  using  SAC2  .  28 

LIST  OF  TABLES 


Table  1.  Triaxial  Loading  under  Drained  Conditions  .  29 

Table  2.  Triaxial  Loading  Under  Undrained  Conditions  .  30 


/ 


v 


Accession  For 

NTIS  GRA&I 
DT lC  TAB 

U;*'u»;:ov>.n 

i  i‘i  cation. 


By - - - 

Distribution/ 


Availability  Codes 
Avail  and/or 
Special 


Dist 


A-i 


INTRODUCTION 


* 


ul 


1 


JjV 

5:;S 

•*/«> 

CL 

•ual 


ill 


l! 


The  objective  of  the  research  reported  herein  was  to  improve  the 
"robustness”  of  the  procedure  used  for  the  numerical  evaluation  of  the 
"bounding  surface”  plasticity  model  for  cohesive  soils. 

The  bounding  surface  model  for  soils  was  developed  and  extensively 
evaluated  under  U.S.  Navy  and  NSF  sponsorship  (1-9).  This  ■  velopment 
included  the  numerical  evaluation  of  the  model  and  the  wri*- of  a  mas¬ 
ter  subroutine  (and  associated  subroutines)  CLAY  that  can  *  incorpo¬ 
rated  into  existing  finite  element  codes  to  supply  the  ncremental 
material  properties  needed  in  the  analysis  of  cohesive  soil  structures 
(2,5,7).  During  the  original  development  phase  emphasis  was  placed  on 
demonstrating  the  capabilities  of  the  model  for  representing  real  soil 
behavior  and  on  the  accuracy  of  the  numerical  evaluation  scheme,  with 
little  attention  given  to  the  robustness  of  the  numerical  evaluation 
(i.e.,  the  robustness  of  CLAY).  CLAY  was  incorporated  into  several  new 
and  existing  finite  element  programs  and  was  used  in  the  analysis  of  a 
number  of  Geotechnical  engineering  problems.  During  the  course  of  these 
studies  several  instances  occurred  in  which  the  numerical  evaluation 
scheme  embedded  in  CLAY  was  far  from  robust. 

In  all  cases  the  problem  involved,  for  one  or  more  elements,  the 
predicted  stress  state  converging  to  a  point  outside  of  the  bounding 
surface.  A  simple  scaling  procedure  had  been  incorporated  into  CLAY  in 
an  attempt  to  avoid  such  an  occurrance,  however,  in  many  cases  it  proved 
to  be  ineffective.  Once  the  stress  state  had  fallen  outside  of  the 
bounding  surface  (at  the  end  of  a  given  solution  ?tep),  in  subsequent 
steps  the  solution  rapidly  deteriorated  and  soon  became  meaningless  and 
very  often  convergence  could  not  be  achieved  (even  to  an  incorrect 
solution),  using  the  original  version  of  CLAY  the  only  remedy  was  to 
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use  smaller  global  solution  increments.  In  order  to  achieve  an  accepta¬ 
ble  solution,  for  many  problems,  the  required  solution  steps  were  exces¬ 
sively  small  and  thus  the  computer  costs  were  excessively  large. 

The  goal  of  the  present  research  was  to  understand  and  to  rectify 
this  problem,  i.e.,  to  develop  a  truely  robust  version  of  CLAY. 
Robustness  requires  that  if  reasonably  sized  solution  steps  are  taken, 
that  accurate  results  be  obtained.  Of  course  in  general,  accurate  one 
step  solutions  cannot  be  achieved  with  a  path  dependent  material  proper¬ 
ties  model.  However,  the  use  of  very  large  steps  leading  to  answers 
which  are  qualitatively  but  not  quantitatively  correct  is  often  useful 
for  a  preliminary  analysis.  Thus,  even  if  ridiculously  large  solution 
steps  are  attempted,  a  robust  routine  should  produce  reasonable  (if  not 
entirely  accurate)  results.  Neither  of  these  conditions  was  met  with 
the  original  version  of  CLAY. 


PRELIMINARIES 

The  bounding  surface  model  can  be  used  to  supply  material  properties 
for  all  solution  schemes  applicable  to  stress  analysis  problems.  The 
most  commonly  used  method  is  the  finite  element  procedure?  of  the  three 
(displacement,  force  and  mixed)  finite  element  formulations  available, 
the  displacement  method  is  most  commonly  used.  Thus,  for  the  remainder 
of  this  report  it  is  assumed  that  the  bouncing  surface  model  is  being 
used  in  conjunction  with  a  displacement  (or  a  reformulated  mixed  [led ) 
finite  element  analysis.  Most  of  the  comments,  however,  also  apply  to 
other  methods,  e.g.,  displacement  formulation  finite  difference 
analyses, etc.  In  the  following  paragraphs  the  finite  element  code,  for 
which  the  material  properties  subroutine  (CLAY)  is  to  supply  material 
properties,  is  called  the  ’’parent”  program. 


The  use  of  a  history  dependent  constitutive  model  (such  as  bounding 
surface  plasticity)  in  a  stress  analysis  program,  in  general,  requires 
some  form  of  incremental  solution  procedure.  In  addition,  if  the  model 
is  nonlinear,  then  in  order  to  be  able  to  employ  reasonably  sized  solu¬ 
tion  steps,  iteration  within  each  step  is  usually  necessary.  The  itera¬ 
tion  of  the  total  solution  by  the  parent  finite  element  program,  will  be 
referred  to  as  ’’global”  iteration. 

For  a  given  point  in  the  body  under  analysis,  and  for  a  given  itera¬ 
tion  K  of  a  given  solution  step  N,  the  role  of  a  ’’properties” 
subroutine,  such  as  CLAY,  is  to  supply  to  the  parent  program  the 
relaionship  between  the  stress  increment  {Aa}  and  the  strain  increment 
{Ae}.  The  bounaing  surface  model  relates  ’’effective”  soil  stress  to 
strain  and,  thus,  throughout  this  section  {a}  and  {Aa}  will  represent 
effective  stresses.  The  fact  that  the  parent  program  maybe  concerned 
with  total  stresses  will  be  addressed  in  a  later  section.  The  collec¬ 
tion  of  points  (locations)  in  the  body  for  which  the  incremental  proper¬ 
ties  are  required  usually  consists  of  all  the  element  centers  or  all  the 
element  integration  points  or  all  the  nodes.  The  following  analysis  is 
in  no  way  dependent  on  which  set  of  points  is  being  used,  however,  for 
simplicity  of  discussion  it  is  assumed  that  a  set  of  properties  is 
required  for  each  element.  The  required  incremental  stress-strain  equa- 
tion  is  usually  written  in  the  form  : 

{Ao}N,K  =  ®N,K-1  *Ae}N,K  +  {A(7o*N,K-1  ^ 

*For  a  ’’force”  (stress)  based  analysis,  the  required  relationship  is 

{Ae}  =  [c]  {Aa}  +  {AeQ} 

This  expression  could  be  found  by  inverting  eq.  1?  however,  it  would  be 
more  economical  to  re-do  the  numerical  evaluation  scheme  in  CLAY  so  as  to 
arrive  directly  at  the  required  expression. 


For  simplicity,  in  the  remainder  of  the  report,  the  above  expression 
will  be  written  without  explicitly  displaying  the  increment  and  itera¬ 
tion  numbers,  i.e. 

{Act}  =  03  {Ae}  +  {Actq}  (2) 

Thus,  the  role  of  the  properties  subroutine,  CLAY,  (for  a  given  glo¬ 
bal  iteration  of  a  given  solution  increment  and  for  a  given  element)  is 
to  provide  the  quantities  03  and  {Aoq}.  For  a  nonlinear  model  they  will 
be  functions  of  {Ao}  and  {Ae},  thus  the  requirement  for  global  itera¬ 
tion.  Typically  the  parent  finite  element  program  supplies  an  estimate 
(from  the  previous  iteration)  of  {Ae}  and  an  estimate  of  {Ao}  is  found 
using  eq.  2  and  the  properties  from  the  previous  iteration  (or  previous 
increment  for  the  first  iteration)?  these  estimates- are  used  by  the  pro¬ 
perties  subroutine  in  the  calculations  of  03  and  {Aoq}-  Now  until  glo¬ 
bal  convergence  occurs,  these  estimates  (for  the  given  element)  of  the 
stress  and  strain  increments  (used  in  the  calculation  of  the  incremental 
properties)  will  not  in  general  satisfy  the  resulting  incremental 
stress-strain  equation,  (2).  Since  this  inconsistency  disappears  as 
global  convergence  takes  place,  it  is  not  absolutely  necessary  to  take 
special  steps  to  rectify  it.  However,  numerical  experimentation  has  demon¬ 
strated  that  overall  computational  advantages  may  be  realized  by  resolv¬ 
ing  it  (5,7).  For  this  purpose,  ’’local  iteration”  can  be  introduced. 
Local  iteration  takes  place  (for  each  element)  within  the  material  prop¬ 
erty  subroutine  (i.e.,  within  CLAY)  and  involves  using  eq.  2  and  the 
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original  strain  estimate  to  calculate  a  new  estimate  of  the  stress* 
which  in  turn  is  used  to  find  a  new  set  of  incremental  properties,  etc. 
The  local  iteration  is  in  addition  to  the  global  iteration.  Since  the 
global  iteration  also  tends  to  remove  the  inconsistency,  the  convergence 
criterion  for  local  iteration  is  typically  less  restrictive  than  the 
global  requirement.  Obviously  the  introduction  of  local  iteration 
increases  the  computational  effort  for  a  given  global  iteration. 
However,  its  use  typically  decreases  the  number  of  global  iterations 
required,  resulting  in  an  overall  reduction  of  computational  effort. 

The  calculations  involved  in  finding  [5]  and  {AaQ}  will  be  discussed 
in  the  next  section. 


MODIFICATIONS  TO  "CLAY” 

The  first  step  was  to  determine  the  cause  of  the  lack  of  robustness 


4 

w 


of  the  numerical  evaluation  scheme  used  in  CLAY.  It  was  found  that 
there  were  two  primary  reasons.  The  first  was  numerical  integration 
error  that  occurred  for  large  (and  sometimes  small)  steps.  The  second 


was  the  inadequacy  of  the  scaling  procedure  used  to  return  a  predicted 

stress  state  to  the  bounding  surface  when  it  fell  outside.  When  the 

solution  steps  were  small  enough  (sometimes  exceedingly  so)  both  of 

these  problems  disappeared  and  the  original  algorithm  functioned  pro 

perly.  The  question  then  was  how  to  improve  the  evaluation  scheme. 

*An  alternative  which  has  not  been  explored  by  the  author  for  the  bound¬ 
ing  surface  model  is  to  maintain  the  initial  stress  estimate  and  itera¬ 
tively  modify  the  strain  estimate.  The  iterative  modification  of  the 
strain  estimate  instead  of  the  stress  estimate  maintains,  at  the  global 
level,  a  compatible  global  displacement  field.  The  rationale  for  this 
choice  is-  the  displacement  continuity  requirement  of  "displacement  for¬ 
mulation"  finite  element  procedures.  The  violation  of  this  condition 
does  not,  however,  necessarily  destroy  the  convergence  characteristics 
of  a  displacement  formulation;  the  effect  it  would  have  in  this  case  has 
not  been  investigated.  It  should  be  noted  that  the  present  local  itera¬ 
tion  scheme  (which  iteratively  modifies  the  stress  estimate)  would  vio¬ 
late  the  equilibrium  requirement  of  a  "force”  formulation  finite  element 
analysis. 
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The  bounding  surface  plasticity  model  leads  to  the  following  rela¬ 
tionship  between  the  stress  and  strain  rates  (see  eq.  4  of  reference 
7): 

{a}  =  Cd]  {e}  (3) 

For  a  precise  definition  of  the  quantity  [o] ,  for  the  bounding  surface 
model,  the  reader  is  referred  to  eq.  5  of  reference  7.  Because  [d]  is  a 
function  of  the  history  of  the  stress  state  over  the  interval,  for  any 
given  solution  step  this  equation  is  nonlinear  (and  thus  its  use,  in 
general,  will  require  iteration}  One  way  of  proceeding  is  to  integrate 
over  the  step  t:yJ  ^->-tN  (in  performing  this  operation  it  is  convenient  to 
think  in  terms  of  time  even  though  the  particular  model  may  be  rate 
independent),  i.e., 


y  {a}  dt  =  y  M{e}  dt 

lN-l  Sl-l 


(4) 


If  it  is  assumed,  for  the  given  solution  step,  that  all  the  strain  com¬ 
ponents  are  proportional  (i.e,  Ae^/Ae^  =  constant  for  t^_1<  t  <  tN), 
then  for  a  rate  independent  model  the  input  history  for  the  interval  can 
be  selected  such  that  all  components  of  {£}  are  constant*  and  are  given 
by  {£}={Ae}/AtN,  hence 


{Aa} 


At  J 
N  t 


Gd 


dt 


N-l 


(5) 


*If  there  are  pore  water  pressure  changes  due  to  water  flow,  then  even 
for  a  rate  independent  model  and  proportionate  strain  increments,  this 
step  will  involve  an  approximation.  This  approximation  occurs  because 
the  assumption  of  a  constant  strain  rate  over  the  interval  would  not 
necessarily  be  consistant  with  the  actual  history  of  the  water  movement. 
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letting 


ra=|-  s  h 

4t”  <-i 

then  gives 

{A cr}  =  [d]  { Ae } 


dt 


(6) 


(7) 


To  this  point  no  approximations  other  than  those  mentioned  above  have 
been  made.  Thus  what  is  needed  is  an  aoeurate  evaluation  of  the  avera- 
age  value  of  [d]  over  the  solution  increment.  Previously,  [2,5,7]  the 
simple  one  step  trapezoidal  rule  was  used,  i.e., 
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where  [d]  n  ^  and  [d]  n  are  the  values  of  [d]  that  correspond  to  the 
stress  and  strain  states  at  the  beginning  and  end  of  the  step,  respec¬ 
tively.  It  must  be  noted  that  in  reality  (]d]n  ^  is  also  a  function  of 
{a}N  (and  hence  {Acr})  because  the  latter  quantity  determines  whether  or 
not  ’’loading”  or  ’’unloading”  occurs  over  the  interval. 

The  previously  cited  lack  of  robustness  has  made  it  evident  that  in 
some  instances  this  simple  numerical  integration  is  not  adequate. 

A  promising  ’’adaptive”  two  point  integration  method  was  tried  but 
failed  to  converge.  Resort  was  then  made  to  the  multi-step  (of  equal 
length)  trapezoidal  rule*. 


*While  the  multi-step  trapezoidal  rule  has  proven  to  be  quite  effective 
it  is  not  the  most  efficient  scheme  available.  Thus,  in  some  future 
research  the  use  of  other  integration  schemes  should  be  explored.  It 
would  seem  to  be  desirable  that  the  end  points  of  the  interval  be 
included  in  the  set  of  integration  points  (because  the  corresponding 
stress  states  explicitly  occur  in  the  incremental  formulation).  Thus 
Simpson’s  rule  and  the  Gaussian  quadrature  method  that  includes  the  end 
points  (Lobatto  quadrature)  would  appear  to  be  the  most  promising  can¬ 
didates. 


The  limit  to  the  accuracy  achieved  with  any  improved  integration 
scheme  is  the  required  assumption  (see  derivation  above)  that  the  strain 
components  vary  proportionately  and  their  rates  are  constant  over  the 
solution  step. 

In  the  following  discussion  the  solution  step  prescribed  by  the  par¬ 
ent  program  (the  finite  element  program  that  is  calling  CLAY)  will  be 
referred  to  as  the  ’’step”.  If  required,  for  accurate  numerical  integra¬ 
tion,  the  algorithm  in  the  new  CLAY  may  further  subdivide  this  step  into 
’’sub-steps”.  These  sub-steps,  however,  will  be  transparent  to  the  par¬ 
ent  program  and  to  the  user  (except  for  increased  computer  cost). 

The  use  of  M  sub-steps  in  the  trapezoidal  rule  leads  to  the  follow¬ 
ing  incremental  properties 

M 

El  =  l  EL  (9) 

_  ,  m 
m=l 

The  incremental  properties  [o]m,  over  sub-step  m,  are  found  from  eq.  8 

where  [d]m  .  and  [d]m,  which  are  the  values  of  [d]  at  times  tM  .  and  t.., 
N-l  N7  N-l  N’ 

are  replaced  by  properties  corresponding  to  times;  (tN  tN_^  +  At^/M) , 
(tN_x  +  AtN/M,  tN_1  +  2AtN/M),  etc.  and  At^  is  replaced  by  At^/M.  The 
strain  at  time  tN  ^  +  m  AtN/M  is  assumed  to  be  (see  previous  paragraphs 
discussing  the  assumptions  made) 


<*>»  *  *  S 


(10) 


The  stress  estimate  at  the  corresponding  time  is  initially  (in  the  first 
local  iteration)  taken  to  be: 


m-1 

{a}m  r  {a}m-l  +  {Aa}m-l  =  {a}N-l  +  {Aa}i  +  <Ao}m-l 


(11) 


That  is,  an  estimate  equal  to  the  value  found  in  the  previous  sub-step 
is  used  for  {Ao>  .  When  local  iteration  is  used,  then  this  is  succes- 
sively  modified  by  replacing  {Ao}^  ^  with  improved  estimates  calculated 
from  eq.  7  (using  {Ae}m  and  [bJm  from  the  previous  local  iteration}. 

The  determination  of  the  number  of  sub-steps  (for  a  given  element) 
to  be  used  in  each  solution  step  is  based  upon  the  following  con¬ 
sideration.  It  was  noted  previously  that  a  symptom  of  inaccurate 
numerical  integration  is  the  prediction  of  a  stress  state  outside  of  the 
bounding  surface.  The  factor  that  measures  this  phenomenon  is  0  which 

relates  the  image  stress  state  a.  .  to  the  actual  stress  state  a.  .  (see 
a  ij  ij 

Figure  1),  i.e. , 

({a}  -  {ac»  =  8({a}  -  {ac»  (12) 

The  image  stress  state  is  defined  by  the  intersection  of  the  projection 

line  with  the  bounding  surface.  The  projection  line  originates  at  the 

projection  center  (see  Figure  1)  and  passes  through  the  stress  state. 

For  the  current  model  the  projection  center  lies  on  the  I  axis  and  is 

located  at  cIQ  (see  Figure  2).  The  size  of  the  bounding  surface  is 

determine  by  the  current  value  of  IQ  (see  [7]).  It  is  seen  from  Figure 

1  and  eq.  12  that  a  value  of  6  of  less  than  1.0  indicates  a  stress  state 

{ct}  which  lies  outside  of  the  bounding  surface. 

In  a  given  global  iteration,  an  initial  attempt  is  made  to  use  one 

step  integration.  If  the  value  of  8,  for  the  calculated  stress  state  at 

the  end  of  the  solution  step  is  greater  or  only  slightly  less  than  1.0 

(>  0.999  ),  then  it  is  assumed  that  no  problem  exists  and  one  step  tra- 

*A  number  of  rather  arbitrary  limits  such  as  this  have  been  used  in  the 
new  version  of  CLAY.  Some  numerical  experimentation  was  done  to  show 
that  the  values  are  adequate,  however,  future  work  should  investigate  a 
range  of  alternatives  in  an  attempt  to  optimize  the  quantities. 


pezoidal  integration  is  used  as  in  the  past  (2,5,7).  When  this  crite¬ 
rion  is  not  satisfied  the  new  CLAY  attempts  to  use  two  sub-steps  in  the 
integration  process.  If  at  the  end  of  either  the  first  or  second  sub¬ 
step  the  3  does  not  satisfy  the  >_  0.999  criterion,  the  suo-step  is  again 
halved  (i.e.  four  sub-steps)  and  the  integration  is  started  again  (using 
one  quarter  of  the  initial  stress  estimate  for  the  beginning  of  the  glo¬ 
bal  iteration).  This  process  is  continued  until  either  the  criterion  is 
met  or  a  limit  of  a  maximum  of  32  sub-steps  is  reached.  In  the  latter 
case  the  solution  process  is  then  continued  as  if  the  limit  on  3  had 
been  satisfied.  Denote  the  number  of  sub-steps  (for  the  given  element) 
arrived  at  by  this  process  as  M  (M  =  1  or  2  or  4  ....32). 

Once,  sub-stepping  is  used,  it  is  used  for  ail  subsequent  global 
iterations  for  that  solution  step  and  for  the  particular  element  in 
question.  As  global  iteration  proceeds,  the  value  of  M  (number  of  sub¬ 
steps)  is  not  permitted  to  decrease  but  it  may  be  increased.*  However, 
in  the  first  global  iteration  of  the  next  solution  step  (for  the  given 
element)  once  again  one  step  integration  is  attemped  (i.e.,  M  is  reset 
to  1).  The  number  of  sub-steps  used  (M),  is  remembered  by  CLAY  (the 
mechanism  will  be  explained  later). 

Within  each  sub-step  local  iteration  is  used  to  determine  the 

value  of  {Aa>.  The  local  iteration  is  continued  until  the  value  of  3 
m 

for  the  end  of  the  sub-step  meets  a  iSG  convergence  criterion  or  a  limit 

initially  a  rather  sophisticated  scheme  was  tried  in  which  the  several 
sub-steps  into  which  the  step  was  divied  could  be  of  varying  length  and 
the  number  of  sub-steps  could  vary  from  one  to  a  set  maximum  as  needed. 
At  the  beginning  of  each  global  iteration  the  number  of  sub-steps  was 
started  at  one,  thus  the  sub-stepping  pattern  could  (and  did)  change 
from  global  iteration  to  global  iteration?  the  resulting  process  had 
poor  convergence  characteristics  (due  to  the  changing  sub-stepping 
pattern)  and  was  thus  abandoned. 

**  Thus  the  control  over  local  iteration  is  no  longer  left  to  the  User 
as  in  the  previous  version  of  CLAY. 


of  a  maximum  of  5  local  iterations  is  reached.  At  least  two  local  inter¬ 
actions  are  always  performed  (the  first  calculation  of  B  is  counted  as 
iteration  number  one).  Because  of  the  local  iteration,  the  initial 
estimate  for  {Ac}  is  not  overly  important  and  a  value  of  zero  could  be 
used  if  desired. 

The  need  for  an  adequate  procedure  for  returning  a  predicted  stress 
state  that  is  outside  of  the  bounding  surface  to  the  surface  still 
exists,  even  though  a  more  accurate  integration  scheme  has  been  adopted. 
Such  a  need  persists  because  of  the  use  of  the  0.999  limit  on  8  instead 
of  1.0,  the  maximum  of  32  placed  on  the  number  of  sub-steps  for  a  given 
solution  step,  and  for  the  intermediate  calculations  in  the  local  itera¬ 
tion  process.  The  limits  on  8  and  M  are  used  to  avoid  excessive  sub¬ 
stepping  and  thus  excessive  computational  cost. 

Classical  ’’radial  return”  (11)  has  been  adopted  to  bring  a  point 
back  to  the  bounding  surface  (while  the  previously  used  scheme  in  CLAY 
was  based  on  the  radial  return  concept,  it  was  only  approximate). 
Whenever  a  stress  state  (at  the  beginning  of  the  step,  or  at  the  end  of 
the  step  or  ore  of  the  sub-steps)  is  found  to  be  outside  of  the  bound  it 
is  scaled  back  to  the  bound  (along  the  line  connecting  it  to  the  projec¬ 
tion  center).  This  scaled  stress  state  is  used  for  the  purpose  of  cal¬ 
culating  the  plastic  modulus. 

The  one  instant  where  the  scaled  value  of  the  stress  is  not  used 
is  in  the  updating  of  the  size  of  the  bounding  surface  as  determined  by 
the  value  of  IQ  (see  eq.  (26)  of  reference  7).  The  extreme  importance 
of  using  the  unsealed  stress  for  this  operation  appears  to  stem  from  the 
fact  that  the  size  of  its  bounding  surface  is  really  controlled  by 
strain  considerations  and  the  strains  are  not  scaled. 


The  scaling  of  the  invariants  of  the  stress  (the  invariants  are 


defined  in  eq.  20  of  (2))  yields i 
1 scaled  =  8(1  -  cIq)  +  cIq 
•^scaled  =  8  J 
Sscaled  = 


(13) 


When  the  stress  state  at  the  beginning  of  the  step  {o}^  ^  is  outside 
of  the  bound,  it  is  of  considerably  more  concern  than  the  corresponding 
problem  for  the  state  at  the  end  of  a  step  or  sub-step.  The  state  at 
the  beginning  of  the  step  is  an  accepted  solution  (to  that  point  in 
time),  whereas  at  the  end  of  a  step  it  is  merely  an  intermediate  predic¬ 
tion  in  the  iteration  process.  The  discrepancy  at  the  beginning  of  the 
step  represents  an  error  which  has  crept  into  the  solution,  whereas  at 
the  end  of  a  step  it  is  only  a  potential  error.  In  an  attempt  to  pre¬ 
vent  error  build-up  (as  incremental  stresses  are  accumulated),  if  the 
stress  state  at  the  beginning  of  the  step  is  found  to  be  oustide-  of  the 
bound,  the  individual  components  are  then  scaled  back  to  the  bound, 


1  •  0  •  j  X 


{a}N-lScaled  =  B{a}N-l  +  Cl  -  3){c}  (14) 

and  a  stress  correlation  vector  {AaQ}  is  calculated 

{Aoq}  =  -  (1-8) ({a}  -  {c})  (15) 

where  {c}  =  j  cIQ  <  1,1, 1,0, 0,0  >T 

This  stress  correction  is  incorporated  into  the  global  analysis  by 
treating  it  as  a  strain  independent  stress  term,  i.e.,  by  modifying  eq. 


The  use  of  the  stress  correction  vector  prevents  a  gradual 
’’straying”  outside  of  the  bounding  surface  for  cases  where  the  state 
should  be  moving  on  the  surface.  Equation  (16)  is  only  used  at  the  glo¬ 
bal  level  because  in  CLAY  the  correction  represented  by  {Aoq}  is 
directly  achieved  by  the  scaling  of  {a}^  i.e.,  all  operations  within 
CLAY  use  eq.  (7). 

For  the  purpose  of  calculating  [□] ,  the  scaling  of  the  stress  states 
back  to  the  bounding  surface  is  of  major  importance,  however,  the  incor¬ 
poration  of  the  stress  correction  vector  {AcrQ}  into  the  global  analysis 
can  be  neglected  with  only  relatively  minor  consequences. 

Situations  can  arise  where  the  one  step  trapezoidal  rule  is  not  ade¬ 
quate,  even  though  it  does  not  result  in  obviously  incorrect  values  of 
8.  An  example  of  such  a  problem  area  is  a  nearly  neutral  loading  situa¬ 
tion  where  inaccurate  integration  may  predict  ’’loading”  when  ’’unloading” 
should  occur  or  vise-versa.  To  avoid  inaccurate  integration,  in 

general,  a  second  sub-stepping  criterion  is  imposed.  The  ratio  |  l"  - 

m ,  n  m  n 
L  /L  is  required  to  be  less  than  0.01,  where  L„  and  L„  are  the  sums  of 

the  absolute  values  of  the  calculated  stress  components  at  the  end  of 
the  increment  with  different  numbers  of  sub-steps  (i.e.  1  and  2  or  2  and 
4,  etc).  This  criterion  will  obviously  always  require  at  least  2  sub¬ 
steps.  As  noted  above,  an  upper  limit  of  32  on  the  number  of  sub-steps 
is  imposed. 

IMPLEMENTATION  OF  CLAY 

Subroutine  CLAY,  along  with  its  supporting  subroutines,  (listed  in 
Appendix  I)  is  intended  to  be  incorporated  into  new  and  existing  finite 
element  (or  finite  difference)  programs  in  order  to  supply  the  incre¬ 
mental  material  properties  for  cohesive  soils  as  predicted  by  the  bound¬ 
ing  surface  plasticity  model.  First  some  general  comments  concerning 


its  use  will  be  made,  followed  by  specific  instructions  for  its  "call”, 
and  finally  some  concluding  remarks  concerning  the  parent  program  and 
the  calculation  of  pore  pressures. 

Subroutine  CLAY  returns  to  the  parent  finite  element  program  tne 
matrices  [5]  and  Od^J  ,  and  the  vectors  {Aa}  and  {Actq}.  The  matrix  [o] 
and  the  vector  { AaQ > ,  as  described  above,  relate  the  incremental  stress 
and  strain  vectors  by  means  of  eq.  2.  The  matrix  [dJ  is  the  tangent 

U 

stiffness  matrix  (12)  and  is  the  value  of  [d]  ,  (eq.  5  of  reference  7) 
at  the  end  of  the  solution  step. 

If  the  global  nonlinear  solution  scheme  uses  successive  substitu 
tion,  only  Gj]  (and  of  course  {A<jq} )  is  needed  (5,i2).  If  the 
Newton-Raphson  method  is  being  used  then  the  Jacobian  can  be  approxi¬ 
mated  by  the  expression  (1-a)  Od]  +  afo^J.  When  a=0  the  procedure 
reverts  to  successive  substitution  and  when  a=l  it  yields  the  ’’tangent 
stiffness”  method*.  For  a  quasi-Newton  method  (13)  neither  fo]  nor  lx>t] 
are  needed  (except  for  initiation  of  the  approximate  Jacobian  and  for 
occasional  updates,  if  used)  and  can  be  ignored.  The  incremental  stress 
{Act},  however,  is  needed  in  the  update  formula  and  is  calculated  in  CLAY 
by  the  previously  described  integration  process  and  is  returned  to  the 
parent  program. 

The  CLAYN  ’’subroutine  package”  consists  of  the  subroutine  CLAY  and 

eight  supporting  subroutines.  A^l  the  subroutines  in  this  package  are 

written  in  FORTRAN  77,  thus  taking  advantage  of  the  structuring  offered 

by  the  language.  In  keeping  with  modern  programming  practice,  the 
£ 

Contrary  to  what  is  suggested  in  (12),  studies  performed  by  the 
authors  for  one  element  problems  have  found  the  successive  substitution 
method  to  show  somewhat  better  convergence  characteristics  than  the  tan¬ 
gent  stiffness  method.  The  successive  substitution  method  can  be  sig¬ 
nificantly  improved  by  using  an  adaptive  convergence  factor  (7). 


design  of  the  package  is  quite  modular  and  facilitates  its  incorporation 
into  new  and/or  existing  finite  element  programs  for  the  analysis  of 
earth  structures. 

To  access  the  CLAY  subroutine  package  a  parent  finite  element  pro¬ 
gram  needs  to  only  call  subroutine  CLAY.  This  is  done  for  each  itera¬ 
tion  of  each  loading  increment  and  for  all  points  in  the  idealized  earth 
structure  where  the  incremental  properties  are  required  (e.g.,  at  ele¬ 
ment  centers  or  at  quadrature  points).  The  call  involves  the  following 
argument  list: 

CALL  CLAY  ( INC ,  ITNO ,  ITYPE ,  KINO ,  LARGE ,  GAMMA.,  PROP ,  STOR ,  SIG ,  EP ,  DSIG ,  DEP ,  U , 
DLTAU , DEAR , DTAN , DS I GO ) 

The  quantities  INC,  ITNO,  ITYPE,  KIND  and  LARGE  are  integer 
variables,  U,  DLTAU,  and  GAMMA  are  real  variables,  and  PROP,  STOR,  SIG, 
EP,  DSIG,  DEP,  DBAR,  DTAN,  and  DSIGO  are  real  arrays  with  the  dimensions 
(21),  (6),  (6),  (6),  (6),  (6),  (6,6),  (6,6)  and  (6),  respectively. 

The  arguments  in  the  call  statement  are  now  discussed  in  detail. 
For  clarity,  the  order  of  discussion  is  changed  slightly  from  that  in 
the  subroutine  call. 

INC;  Global  solution  increment  number  (the  first  increment  must  be 
numbered  1). 

ITNO:  Global  iteration  number  (the  first  iteration  in  each  increment 

must  be  numbered  1). 

ITYPE:  A  flag  indicating  which  form  of  the  bounding  surface  is  to  be 

used  in  the  analysis  (further  details  regarding  the  possible 
forms  of  the  surface  are  given  in  Appendix  III). 

PROP:  A  one-dimensional  array  containing  the  values  of  the  model 

parameters.  At  the  point  in  the  idealized  earth  structure  for 
which  the  incremental  material  properties  are  sought,  these 


parameters  characterize  the  bounding  surface  plasticity  model 
for  the  soil.  The  parent  finite  element  program  must  read  and 
store  the  values  of  the  model  parameters  for  each  different 
type  of  soil  in  the  structure  being  analyzed,  and  then,  for 
each  call  to  CLAY,  present  the  appropriate  values  for  the  soil 
at  the  point  in  question. 

The  model  parameters  are  stored  in  the  PROP  array  in  the 
following  order*:  A,  <,  Mc,  Mg/Mc ,  v  or  G,  r,  ,-  Pgtm,  Rc, 

V  T»  W  Ae/Ac’  c’  s’  m’  hc’  W  h2’  a’  That  is’ 
PR0P(1)  =  A,  PR0P(21)  =  w,  etc.  After  the  model  parameters  are 

read  into  the  parent  program,  but  before  the  first  call  to 

CLAY,  the  values  for  M  and  p.  must  be  converted  to  associated 

c  x# 

quantities  in  invariant  stress  space.  This  is  achieved  by  mul¬ 
tiplying  these  parameters  by  1//27  and  3,  respectively.  It  is 
suggested  that  subroutines  RPROP  and  TCHECK,  listed  in  Appendix 
II,  be  incorporated  into  the  parent  finite  element  program  for 
the  purpose  of  reading,  echo  printing  and  scaling  the  material 
parameters.  The  input  instructions  for  RPROP  are  given  in 
Appendix  III  along  with  a  brief  discussion  of  the  new  single 
surface  option  for  the  bounding  surface. 

STORi  This  one-dimensional  array  is  used  to  store  the  values  of  cer¬ 
tain  quantities  (e.g.,  internal  variables)  which  describe  the 
current  state  of  the  soil  and  parameters  related  to  the  evalua¬ 
tion  of  the  model  such  as  the  number  of  sub-steps  M.  For  a 
given  step  in  the  analysis  these  values  are  unique  to  the  point 
in  the  structure  under  consideration  (e.g.,  element  centers). 
At  the  beginning  of  the  analysis  the  parent  program  must  ini¬ 
tialize,  for  each  point  in  question,  the  values  of  ST0R(1)  and 
*A  detailed  discussion  of  the  model  parameters  is  given  in  ref.  6,7. 


ST0R(5).  STOR(l)  must  contain  the  initial  value  of  the  effec¬ 
tive  preconsolidation  pressure  (IQ  =  3pQ),  while  the  initial 

value  of  the  total  void  ratio  (eQ)  must  be  stored  in 
* 

5T0R(5J.  These  values  are  not  read  by  the  material  parameter 
input  subroutine  (RPROP)  because  they  will,  in  general,  vary 
from  point  to  point  within  the  deposit,  even  for  a  homogeneous 
soil.  After  eacn  call  to  CLAY,  the  values  in  5T0R  must  be 
stored  by  the  parent  program  for  each  point  in  the  earth  struc¬ 
ture  for  which  the  incremental  properties  are  required.  Prior 
to  each  call  to  CLAY  the  appropriate  values  for  the  point  in 
question  must  be  retrieved  from  storage  (i.e.,  from  a  two- 
dimensional  array  or  from  a  disk  file)  and  then  presented  to 
the  CLAY  subroutine  through  the  CALL. 

SIGi  This  one-dimensional  array  contains  the  components  of  the 

stress  at  the  beginning  of  the  increment  and  must  be  supplied 

to  CLAY  by  the  parent  program.  It  corresponds  to  the  vector 
{a}^  Compressive  normal  stresses  are  considered  to  be 

positive. 

EP:  The  components  of  the  strains  at  the  beginning  of  the  increment 

are  stored  in  this  one-dimensional  array  and  must  be  supplied 
to  CLAY  by  the  parent  program.  It  corresponds  to  the  vector 
{e>N  Compressive  normal  strains  are  considered  to  be  posi¬ 
tive. 

OSIGi  This  one-dimensional  array  contains  an  estimate  of  the 

stress  increment.  It  corresponds  to  the  vector  {Ao}^  and  is 
calculated  and  returned  to  the  parent  program  by  CLAY. 

& 

Note:  When  the  new  version  of  CLAY  is  introduced  into  finite  element 
programs  SAC2  and  SAC3  (9),  the  storage  of  the  initial  void  ratio  in 
ST0R(7)  must  be  changed  to  ST0R(5)  in  Subroutines  GEOM  and  STIFNS. 

**Whether  these  are  to  be  ’’total”  or  ’’effective”  stresses  is  discussed 
later. 
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The  estimate  of  the  strain  increment  is  contained  in  this  one¬ 


dimensional  array.  It  corresponds  to  the  vector  {Ae)N  and  is 
supplied  by  the  parent  program. 

|J:  This  variable  represents  the  core  pressure  at  the  beginning  of 

the  increment  N-l.  It  corresponds  to  the  quantity  uM  * 

DLTAUi  The  estimate  of  the  pore  pressure  increment  is  represented  by 
this  variable.  It  corresponds  to  the  quantity  Au  *. 

GAMMA:  This  variable  is  the  combined  bulk  modulus  (r)  for  the  soil 

particles  and  pore  fluid  (14)*. 

KIND:  The  value  assigned  to  this  flag  is  determined  by  how  saturated 

conditions  are  modelled  in  the  parent  program. 

A  value  of  zero  is  required  when  the  special  (mixed)  formu¬ 
lation  for  incompressible  and  nearly-incompressible  solids 
(10,15)  is  used.  In  such  instances,  the  pore  pressure  is 
treated  as  a  primary  dependent  variable  at  the  global  level. 
Thus,  the  parent  program  calculates  uN  and  Au^  1  at  the  glo¬ 
bal  level,  whether  or  not  the  parent  program  is  required  to 
send  them  to  CLAY  (in  u  and  DLATU)  depends  on  whether  it  sends 
total  or  effective  stresses  in  SIG  and  DSIG  (these  possibili¬ 
ties  are  explored  in  greater  detail  later).  For  such  a  for¬ 
mulation,  the  value  of  r  (GAMMA)  will  be  required  by  the  parent 
program  and  is  supplied  by  CLAY. 

When  a  conventional  displacement  formulation  for  compressi¬ 
ble  solids  is  used  for  idealized  undrained  conditions  (14),  (or 
for  unsaturated  conditions  where  r  is  equal  to  zero),  KIND  is 
set  equal  to  one.  In  such  cases,  the  only  primary  global  vari¬ 
ables  are  the  displacements;  the  pore  pressure  is  treated  as  a 

* 

Whether  this  quantity  is  supplied  by  the  parent  program  or  by  CLAY  is 

discussed  in  the  explanation  of  the  variable  KIND. 


secondary  dependent  variable.  The  CLAY  subroutine  calculates 


the  values  of  uN  ^  (U)  and  Au^  (DLTAU),  and  supplies  them  to 
the  main  program  for  printing. 

OBAR  and 

0TAN»  These  arrays  contain  the  estimates  of  the  incremental  stress- 


strain  properties 


(see  eq.  2)  and  the  ’’tangent  stiffness 


matrix”  [d^J  respectively.  For  KIND  =  1  these  matrices  have 
been  augmented  by  adding  r  to  the  matrix  elements  in  the  upper 
left  3x3  corner  (more  details  are  given  later;  see  also 
7,14). 

Subroutine  CLAY  calculates  [d]  and  [d J  and  returns  them 
to  the  parent  program.  It  is  desirable,  but  not  absolutely 
necessary,  that  the  parent  program  ’’remembers”  the  last  [o] 
calculated  for  the  point  in  question  and  returns  it  to  CLAY  in 
the  next  call.  This  is  similar  to  the  requirement  for  STOR. 
For  the  first  iteration  of  the  first  solution  step  (and  at 


other  instances  for  which  the  last  calculated 
available)  it  should  be  set  to  zero. 


is  not 


DSIGOi  This  array  contains  the  stress  correction  vector  of  eq.  2,  i.e. 
{Aoq}. 

LARGE «  The  value  of  this  flag  depends  on  the  definitions  used  for  the 
stresses  and  strains  in  the  analysis.  If  engineering  measures 
of  stresses  and  strains  are  used,  LARGE  is  set  equal  to  zero. 
If  true  stresses  and  logarithmic  (natural)  strains  are  used, 
LARGE  is  set  equal  to  one. 

Although  the  CLAY  subroutine  package  determines  three-dimensional 
incremental  stress-strain  properties,  it  can  also  be  used  to  supply  pro¬ 
perties  for  two-dimensional  analyses.  For  a  three-dimensional  analysis 
the  ordering  of  the  components  in  the  stress  (and  strain)  vector  is 


In  order  to  use  CLAY  in  conjunction  with  a  two-dimensional  analysis 
the  strain  and  stress  vectors  in  the  parent  program  must  be  expandea 
into  this  form.  For  example,  for  an  axisymmetric  analysis  the  stress 
and  strain  components  are  ordered  in  the  following  manner* 

<ar,ae’az’Tre*  °’°’>T  and  <er,£8’£z,Yr0’  °*  0>T 

The  indicated  zero  values  must  be  supplied  by  the  parent  program  for  the 
vectors  {a}^,  {e}^,  and  ^e^N*  The  ^ncrementa^  stress-strain  proper¬ 
ties  of  interest  are  located  in  the  upper  left  4x4  corners  of  the  6x6 
DBAR  and  DTAN  arrays  returned  by  the  CLAY  subroutine,  etc. 

For  plane  strain  conditions  the  stress  and  strain  vectors  are 
ordered  in  the  following  manner* 

«V  °y»  <i2>  Txy>  0,  0>T  and  <ex,  cy,  0,  yxy,  0,  0>T 

The  indicated  zero  values  must  be  supplied  by  the  parent  program  for  the 
appropriate  arrays.  It  is  seen  that  the  parent  program  must  calculate 
the  stress  normal  to  the  plane  of  the  body  (i.e.,  a  ). 

As  described  to  this  point,  subroutine  CLAY  provides  an  incremental 
relationship  between  the  strain  increment  •  and  the  effective  stress 
increment.  (Of  course,  for  ideal  drained  conditions  total  and  effective 
stresses  are  identical).  If  CLAY  is  to  be  incorporated  into  a  parent 
finite  element  program  which  only  requires  this  relationship,  then  no 
further  considerations  are  needed.  For  such  an  application  the  "call 
parameters”  KINO,  U  and  DLTAU  are  assigned  values  of  0,  0.0  and  0.0 
respectively,  and  the  material  property  GAMMA  is  ignored.  The  arrays 
SIG  and  DSIG  are  then  effective  stresses. 


Saturated  soils  under  ideal  undrained  conditions  (no  movement  of 
pore  fluid)  behave  as  nearly  incompressible  solids  and  are  often  approx¬ 


imated  as  perfectly  incompressible  solids.  For  a  perfectly 
incomDressibie  solia  a  conventional  displacement  analysis  can  not  be 
used  (10).  The  most  popular  means  for  handling  this  problem  are  the 
use  of  a  ’’reformulated”  analysis  (i.e.  a  ’’mixea”  approach)  (10)  ana  tne 
use  of  a  penalty  approximation  of  the  incompressibility  (16).  The 
reformulated  analysis  only  requires  a  relationship  between  effective 
stress  and  strain  and  hence  the  comments  of  the  previous  paragraph 
apply.  The  penatly  method  in  reality  considers  the  material  to  be  only 
nearly  incompressible  and  is  dealt  with  in  the  next  paragraoh. 

It  is  often  assumed  that  the  slight  compressibility  of  the  pore 
fluid  and  the  soil  particles  can  be  modeled  by  the  simple  linear  rela¬ 
tion: 

AVAq  =  -  TAu 

Where  AVAq  is  the  incremental  change  in  volume,  Pis  a  constant  bulk 
* 

modulus  for  the  pore  fluid-soil  particles  and  Au  is  the  change  in  pore 
pressure.  (The  approximation  of  incompressible  behavior  by  means  of  a 
.^-'''penalty  formulation  introduces  this  same  equation  where  r  is  interpreted 
as  the  penalty  number  (16)).  Typical  Twill  be  very  large  compared  to 
the  deviatoric  stiffness  and  thus  the  soil  will  behave  as  a  ’’nearly 
incompressible  solid.”  The  analysis  of  nearly  incompressible  solids  is 
notoriously  difficult  (10,17)  and  requires  some  special  consideration. 
The  most  commonly  used  schemes  are  the  conventional  displacement  formu¬ 
lation  with  ’’selected-reduced”  integration  (16)  and  the  previously  men¬ 
tioned  reformulated  (mixed)  analysis  (10). 

* 

The  possibility  of  treating  T  as  a  variable  in  an  attempt  to  model  par¬ 
tially  saturated  conditions  is  beyond  the  scope  of  this  study. 
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When  a  reformulated  analysis  is  used  for  a  nearly  incompressible 
material  the  parent  program  needs,  in  addition  to  the  effective  stress- 
strain  relationship,  the  volume  change-pore  pressure  relationship.  For 
this  latter  purpose  CLAY  supplies  the  value  of  r  (GAMMA  of  the  call)  to 
the  main  program.  In  a  reformulated  analysis  the  pore  water  pressure  is 
directly  calculated  in  the  global  finite  element  analysis.  If  it  is 
desirea  that  the  parent  proqram  send  and  receive  total  stress  informa¬ 
tion  (SIG  and  DSIG)  to  CLAY  then  it  must  also  supply  (in  the  call)  the 
appropriate  pore  pressure  values  (U  and  DLTAU)  to  CLAY.  If,  instead, 
effective  stresses  are  used  then  u  and  DLTAU  are  supplied  as  zero. 

For  a  conventional  displacement  analysis  of  idealized  undrained  con¬ 
ditions  where  the  water-soil  particle  compressibility  is  included,  the 
material  stiffness  is  augmented  by  the  volume  compressibility  due  to  r, 
i.e.,  r  is  added  by  CLAY  to  each  member  of  the  upper  left  3x3  corner 
of  GO  and  (jD^J  (OBAR  and  OTAN  of  the  call).  For  such  conditions  the 
call  parameter  KIND  is  set  equal  to  1,  and  SIG  and  DSIG  are  total  stress 
vectors.  Subroutine  CLAY  calculates  the  resulting  pore  pressure  at  the 
beginning  of  tne  interval  (U)  and  the  incremental  change  (DLTAU)  and 
returns  them  to  the  main  program  for  printing  purposes. 

Either  of  the  two  methods  discussed  for  undrained  conditions  are 
valid  for  ideal  drained  conditions  if  r  is  set  equal  to  zero. 

The  most  important  class  of  problems  is  when  the  soil  is  saturated 
and  there  is  time  for  flow  of  the  pore  fluid  to  take  palce  (due  to  non- 
homogenous  stress  conditions).  Such  situations  require  the  solution  of 
the  coupled  flow-stress  analysis  problem  (a  number  of  references  are 
given  in  (8)).  Such  an  analysis  requires  the  incremental  strain- 
effective  stress  relationship  which  can  be  supplied  by  CLAY.  The. param¬ 
eter  KIND  is  set  equal  to  0.  The  stress  vectors  (SIG  and  DSIG)  are 


either  effective  or  total  depending  on  whether  u  and  DLTAU  are  supplied 
to  CLAY  as  zeros  or  as  the  actual  pore  water  pressure  values  (as  calcu¬ 
lated  by  the  parent  program  from  the  coupled  analysis). 


The  subroutines  making  up  the  CLAY  subroutine  package  are  briefly 

described  below.  Listing  of  the  subroutines  are  given  in  Aopendix  I. 

CLAY:  This  is  the  main  subroutine  of  the  package.  The  parent  finite 

element  program  must  call  this  subroutine  in  the  manner  dis¬ 
cussed  in  the  previous  section.  The  substepping  integration 
strategy  is  determined  and  the  local  iteration  is  conducted  in 
the  subroutine. 

BOUND 1 :  In  this  subroutine  the  quantities  8,  F,  F,j,  F,a  and  Kp  are 
computed  for  the  form  of  the  bounding  surface  consisting  of  a 
single  ellipse  (further  details  regarding  this  new  form  of  the 
surface  are  given  in  Appendix  III  and  in  (18));  i.e.,  for 
ITYPE=1  (refer  to  the  discussion  of  the  call  to  CLAY). 

B0UND3 :  Values  for  the  quantities  listed  above  are  computed  for  the 
form  of  the  bounding  surface  consisting  of  two  ellipses  and  a 
hyperbola  (7);  i.e.,  for  ITYPE=3. 

ELASTC:  The  elastic  contribution  to  the  [b]  array  is  computed  in  this 
subroutine. 

GETH:  The  value  of  the  hardening  function  is  determined  in  this 

subroutine.  Further  details  regarding  forms  of  the  hardening 
function  are  given  in  Appendix  III. 

INVAR:  The  values  of  the  three  stress  invariants  used  in  the  formula¬ 

tion  are  computed  in  this  subroutine. 

LODFUN:  The  values  of  the  plastic  modulus,  K  ,  and  of  the  loading 
index,  L  (eq.  (7)  in  reference  (7)),  are  computed  in  this  sub¬ 
routine. 
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PLASTCj  The  plastic  contribution  to  the  M  array  (eq.  (5)  in  reference 
(7))  is  computed. 

SKALEs  In  this  subroutine  if  the  stress  state  lies  outside  the  bound¬ 
ing  surface,  then  using  eqs.  (13),  it  is  scaled  back  to  it. 

EXAMPLES 

A  number  of  examples  have  been  analyzed  in  order  to  demonstrate  the 
improvements  in  ’’robustness”  of  CLAY.  The  first  series  involved  incor¬ 
porating  the  new  CLAY  into  program  EVAL  (7)  and  solving  single  element 
examples.  For  those  cases  in  which  CLAY  had  previously  exhibited  poor 
behavior  the  improvement  was  dramatic.  In  Tables  1  and  2  comparisons 
are  made  of  the  results  from  two  of  these  analyses  for  varying  numbers 
of  global  steps  to  reach  a  given  final  loading  condition,  while  the 
results  obtained  with  the  new  CLAY  for  one  very  large  global  step  are 
not  entirely  accurate,  they  are  not  unreasonable  as  was  the  case  with 
the  old  version  of  CLAY.  A  careful  study  of  the  predicted  pore  water 
pressures  (u)  in  Table  2  for  the  new  CLAY  reveals,  an  interesting  pheno¬ 
menon,  i.e.,  the  convergence  with  increasing  numbers  of  global  incre¬ 
ments  is  not  entirely  smooth.  This  behavior  can  be  observed  by 
comparing  the  results  with  2,  4  and  8  global  steps.  This  phenomenon  is 
caused  by  the  adaptive  sub-stepping  scheme  used  in  the  new  CLAY  which 
can  result  in  a  greater  total  number  of  sub-steps  even  though  the  number 
of  global  steps  is  less. 

The  second  series  involved  using  the  old  and  new  versions  of  CLAY  in 
the  plane  strain  finite  elementr  program  SAC2  and  analyzing  two  ideal¬ 
ized  footing  problems.  Properties  for  a  typical  unsaturated  clay  were 
used  for  the  soil.  Results  from  one  of  these  examples  are  given  in 
Figure  3.  The  improvement  in  predictions  for  large  step  sizes  is 
clearly  evident. 
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CONCLUSIONS 

The  numerical  scheme  used  to  evaluate  the  bounding  surface  model  for 
cohesive  soils  has  been  successfully  modified  in  order  to  substantially 
improve  its  robustness.  The  modifications  include  the  introduction  of 
sub-stepping  when  necessary  and  a  more  accurate  radial  return  algorithm. 

It  was  demonstrated  by  a  number  of  examples  that  the  new  algorithm 
is  quite  robust.  That  is  for  reasonably  sized  solution  steps  it  gives 
accurate  results  and  for  very  large  steps  it  is  convergent  and  gives 
reasonable  accuracy.  The  one  question  that  remains  to  be  explored  is 
the  impact  of  these  modifications  on  the  computational  cost.  The  cost 
of  the  material  model  evaluation  (as  compared  to  other  costs  such  as 
equation  solution)  for  large  finite  element  analyses  of  earth  structures 
must  be  determined. 
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Table  1.  Triaxlal  Loading  under  Drained  Conditions 
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ooo  ooo  oooo  o  o  n 


SUBROUTINE  CLAY  ( INC , ITNO , ITYPE , KIND , LARGE , GAMMA , PROP , STOR , SIGTM , 

$  SPM , DSIGM , DEPM , UB.DLTAU, DEAR , DTAN , DSIGO ) 

'■  Subroutine  to  evaluate  Yannis  Daf alias'  bounding  * 

*  surface  plasticity  model  for  clay  soils.  * 

*  Fortran  77  version,  Prepared  by  L.R.  Herrmann  and  V.Kaliakin  * 

*  at  the  University  of  California,  Davis  Campus.  * 

»  • 

*  "Robust”  version  using  radial  return  and  3ub-increaenting  -  11/85  0 

INTEGER  I , J , K , IT , INC , ITNO , ITYPE , KIND , LARGE , LOIT , LOITMX ,11(6) 

REAL  PROP( 1 ) ,ST0R(6 ) , SIGTM(6 ) ,EPM(6 ) ,DSIGM(6 ) ,DEPM(6 ) , DBAR(6 ,6 ) , 

$  DTAN (6,6), SIGB ( 6 ) , EPB ( 6 ) , DSIG ( 6 ) , DEP ( 6 ) , DEPT ( 3 , 3 ) >  SB ( 3 , 3 ) > 

$  SE(3, 3 ) ,DLTA(3 ,3 ) ,DB(6 , 6 ) ,DSIGQ(6 ) ,SIGEM(6 ) ,SIGBP(6 ) , 

$  UB , DLTAU , GAMMA , SMALL , DIL , DDIL , VOIDB , VOIDE , XIB , XIE , XJB , X JE , 

$  SCUBEB , SCUBEE , SIN3AB , SIH3AE , COS3AB , COS3AE , BULKB , 3ULKE , GB , GE , 

$  XI OB , XIOE , XIL , XIBSV , XIESV , BETA , ER , REF , CONV , PRT , RTS , 

$  TEMP 1 , TEMP2 , TEMP 3 , TEMP4 

C 

DATA  11/11, 22, 33, 12, 13, 23/,  DLTA/1 .0, 3*0,0, 1 .0, 3*0.0, 1 .0/ 


A  number  of  arbitrary  limits  are  defined 

SMALL=0.0Q01 *PR0P(8 )  !  small  value  for  stress  invariant 
BETALM=0.999  !  limit  on  "beta"  (how  far  can  be  out  of  bound) 

STPMIN=1 ,/32.  +  .01  !  limit  of  32  sub-steps 

L0ITMX=5  !  limit  on  number  of  local  iterations 

C0NV=0. 01  !  convergence  limit  for  local  iterations 

and  sub-stepping 


Initialize  history  for  first  increment 

IF (INC  .EQ.  1  .AND.  ITNO  .EQ.  1)  THEN 

ST0R(2)=ST0R( 1 )  !  size  of  bounding  surface  Io 

ST0R(3)=UB  !  pore  water  pressure(for  non-ref ormulatea ) 

ST0R(6)=1.Q  !  size  of  sub-steps 

END  IF 

Update  history  for  new  increment 

IF (INC  .GT.  1  .AND.  ITNO  .EQ.  1)  THEN 

ST0R( 1 )=ST0R(2)  !  size. of  bounding  surface  Io 

STOR(3)=STQR(4 )  !  pore  water  pressure(for  non-ref ormulated) 

ST0R(6)=1.0  !  size  of  sub-steps 

END  IF 

Convert  from  total  stress  formulation  to  effective  stress 
GAMMA=PROP(6 )  l  bulk  modulus  for  water  &  soil  particles 


IF(XIND  .  EQ.  1)  THEN  !  unreformula ted  analysis 
UB=STOR(3 ) 

DLTAU=GAMMA*(DEPM(  1  )+DEPM(2)-rDEPM(5 ) )  ichange  in  pore  pressure 
ST0R(4)=UB  +  DLTAU  !  store  pore  press,  (unreforaulated  analys) 

DO  5  1=1,3  1  remove  stiffness  due  to  gamma 

DO  4  J=1 ,3 

DBAR ( I , J ) =DBAR ( I , J )  -  GAMMA 

4  CONTINUE 

5  CONTINUE 
ENDIF 

DO  10  1=1,3  !  convert  to  effective  stresses 

SIGEM(I )=SIGTM(I )  -  U3 

SIGBP(I)=SIGEM(I)  !  initial  guess  of  stress  for  end 

SIGEM(I+3 )=SIGTM(I+3 )  !  of  increment 

3IGBP(Ii-3  )=SIGEM(I+3 ) 

10  CONTINUE 

Calculate  the  initial  estimate  for  the  incremental  stress 

DO  30  1=1,6 
TEMP=0. 0 
DO  20  J=1  ,6 

TEMP=TEMP  -r  DBAR  ( I ,  J )  *DEPM  ( J ) 

20  CONTINUE 

DSIGM(I)=TEMP 
30  CONTINUE 

RTS=ST0R(6)  !  recall  sub-step  size  used  in  previous  iteration 

40  CONTINUE  !  return  point  when  change  sub-step  size 

Transfer  stress  , strain  and  pore  water  pressure  to  local  arrays  so 
as  not  to  disturb  the  values  brought  in  from  the  parent  program 

DO  100  1=1,6 

SIGB ( I ) =SIGEM ( I ) 

DSIG ( I ; =DSIGM ( I ) *RTS 
SPB(I)=£PM(I) 

DEP(I)=DEPM(I)*RTS 
DSIG0(I)=0.0 
DO  50  J=1 ,6 

DBAR(I , J)=0. 0 
DTANCI. J)=0.0 
50  CONTINUE 
100  CONTINUE 

Determine  3-dimensional  incremental  properties. 

FACT0R=0.0  I  Tangent  properties  at  start  of  increment 

XI0B=ST0R( 1 ) 

CALL  INVAR  ( FACTOR, STOR, SIGB, DSIG, EPB.DEP, DEPT, V0IDB, LARGE, 
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$  SMALL , XIB , XJB , DIL , DDIL , SB , SCUBEB , SIN3AB , C0S3AB , 

0  DLTA , II } 

XIBSV=XIB  !save  unsealed  value  for  finding  change  in  lo 
CALL  ELASTC  ( PROP, VOIDB, XIB, DB,3ULKB,GB) 

CALL  PLA3TC  ( ITYPE , PROP , DEPT , VOIDB , XIB , XJB , XIOB , DDIL , SB , SCUBEB , 

$  SIN3AB , C0S3AB , DB , BULXB , GB , DLTA , II , INC , ITNO , BETA ) 

C 

C  If  stress  at  start  of  increment  outside  of  bound  move  it  oacic 
C  and  calculate  stress  correction  vector 

C 

IF (BETA  .LT.  1.0)  THEN 
DO  200  1=1,3 

DSIG0(I)=-(1 .0  -  BETA)«(SIGB(I)  -  XIB*PROP( 14 )/3 . 0 ) 

SIGB(I)  =SIGB(I)  +  DSIGO(I) 

DSIG0(I+3)=-(1 .0  -  BETA )*SIGB( 1+3) 

3IGB( 1+3 )  =SIGB(I+3)  +  DSIG0(I+3) 

200  CONTINUE 
END  IF 

Sub-incrementing  (if  needed)  integration  loop  over  solution  step 

PRT=0 . 0 
C 

300  CONTINUE  !  Sub-incrementing  loop 

C 

BETAL=0. 0  I  initialize  memory  for  beta 

DO  325  L0IT=1 ,L0ITMX  !  Local  iteration  on  stress  sub-increment 
C 

FACT0R= 1 . 0  !  Tangent  properties  at  end  of  sub-increment 

CALL  INVAR  ( FACTOR, ST0R,SIG3,DSIG,EP3,DEP, DEPT, VOIDE, LARGE, 

$  SMALL , XIE , XJE , DIL , DDIL , SE , SCUBEE , SIN3AE , C0S3AE , 

$  DLTA, II) 

XIESV=XIE  J  save  unsealed  value  for  finding  change  in  Io 

C 

CALL  ELASTC  (PROP, VOIDE, XIE, DTAN , BULXE , GE ) 

C 

XIL=PR0P(7)  !  calculate  size  of  bounding  surface 

TEMPI = 1 . 0/ ( PR0P( 1 )  -  ?R0P(2 ) ) 

TEMP2= (XIESV  -  XIBSV)/3.0 
IF (XIOB  .GE.  XIL)  THEN 

XI0E=XIQB*EXP(TEMP1 *0 . 5*( (VOIDB  +  VOIDE)*DDIL 
$  -  (VOIDB/BULXB  +  V0IDE/BULKE)*TEMP2) ) 

ELSE 

XIQE=XI0B  +  TEMPI *XIL*Q.  5#(  (V0IDE+V0IDB)*DDIL 
$  -  (VOIDB/BULXB  +  V0IDE/BULXE)»TEMP2) 

END  IF 

CALL  PLASTC  (ITYPE, PR0P,DEPT, VOIDE, XIE, XJE, XI0E, DDIL, SE, SCUBEE, 

$  SIN3AE , C0S3AE , DTAN , BULXE , GE , DLTA , II , INC , ITNO , BETA ) 

C 

DO  320  1=1,6  !  estimate  of  stress  change  over  sub-increment 

TEMP=0.0 


ooooo  n  ooooo 


TEMP=TEMP  +  0.5*(DB(I,J)  +  DTAN(I, J) )*DEP(J) 

310  CONTINUE 

DSIG(I)=TEMP 

320  CONTINUE  !  check  convergence 

IF ( ABS ( 3ETA-BETAL ) /MAX ( BETAL , BETA )  .LT.  CONV)  GOTO  326 
BETAL=BETA  I  update  memory  of  beta 

325  CONTINUE 

326  CONTINUE 


Check  to  see  if  sub-incrementing  is  sufficiently  fine 

IF (BETA  .LT.  BETALM  .AND.  RTS  .GT.  STPMIN)  THEN 

RTS=RTS*0.5  !  halve  the  interval 

GO  TO  40  !  start  over  on  the  step  with  smaller  sub-steps 

END  IF 

PRT=PRT  +  RTS  !  update  at  end  of  sub-step 

BULKB=BULKE 

XI3SV=XIESV 

XIB=XIE 

VOIDB=VOIDE 

XIOB=XIOE 

DO  400  1=1,6 

SIGB ( I ) =SIGB ( I )  +  DSIG(I) 

EPB(I)  =EPB(I)  +  DEP(I) 

DO  390  J=1 ,6 

DBAR ( I , J ) =DBAR ( I , J )  +  0.5*(DB(I,J)  +  DTAN(I, J) )*RTS 
DB ( I , J ) =DTAN ( I , J ) 

390  CONTINUE 
400  CONTINUE 

IFCPRT  .LT.  0.99)  GOTO  300  !  solution  increment  not  yet  complete 

Check  to  see  if  sub-stepping  is  sufficiently  fine  to  give 
accurate  stress  predictions 


ER=Q. 0 
REF=0. 0 
DO  450  1=1,6 

ER=ER  +  ABS(SIGB(I)  -  SIGBP(I)) 

REF = REF  +  ABS(SIGB(I) ) 

SIGBP(I)=SIGB(I)  !  remember  stress  with  previous 

450  CONTINUE  !  step  size 

IF(ER/REF  .GT.  CONV  .AND.  RTS  .GT.  STPMIN)  THEN 
RTS=RTS#0.5 

GOTO  40  I  start  over  on  the  step,  only  with 

END  IF  I  smaller  sub-steps 

IF (BETA  .LT.  1 . 0 )WRITE ( 4 , * ) ' ITN0= ' , ITNO , ' BETA= ' , BETA 


C 


ST0R(2 : =XI0E 


!  store  size  of  bounding  surface 


n  o  o  o  o  o 


ST0R(6)=RTS  !  store  sub-step  size 

Calculate  incremental  stress  and  return  it  to  the  main  program 

DO  500  1=1,6  !  stress  increment 

DSIGM ( I ) =SIGB ( I )  -  SIGEM(I) 

500  CONTINUE 

Conversion  from  effective  stress  to  total  stress 

IF(KIND  .EQ.  1}  THEN  !  for  non-reformulated  analysis 

DO  560  1=1,3  l  add  in  stiffness  due  to  gamma 

DO  555  J=1 ,3 

DBAR ( I , J ) =DBAR ( I , J )  +  GAMMA 
DTAN(I, J)=DTAN(I, J)  +  GAMMA 
555  CONTINUE 

560  CONTINUE 
END  IF 

DO  600  1=1,3  !  add  pore  pressure  to  effective  stress 

DSIGM ( I ) = DSIGM ( I ) +DLTAU 
600  CONTINUE 
RETURN 
END 
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SUBROUTINE  INVAR  ( FACTOR, STOR,SIG,DSIG,EPB,DEP, DEPT, VOID, LARGE, 
$  SMALL , XI , X J , OIL , DDIL , S , SCUBE , SINS A , CQS3 A , 

$  DLTAjII ) 

Subroutine  to  compute  values  ox"'  invariants 

INTEGER  I , J , K , N , IK , LARGE ,11(1  ) 

REAL  STOR(6),SIG(6),DSIG(6),£PB(6),DEP(6),DEPT(3,3), 

$  S(3, 3), DLTA(3, 3), VOID, SMALL, XI, XJ.DIL, DDIL, SCUBE, 3IN3A, 

$  C0S3A, ARB, FACTOR, TEMPI 

PARAMETER  (yes=1 ) 

DATA  ARB/ 1000.0/ 

DIL=0.0 
DDIL=0.0 
DO  20  1=1,3 

DIL  =DIL  +  EPB(I) 

DDIL=DDIL  +  DEP(I) 

20  CONTINUE 
XI=0.0 
XJ=0.0 
SCUBE=0.0 
SIN3A=0.0 

Calculate  first  stress  invariant 
DO  100  1=1,3 

XI=XI  +  SIG(I)  +  FACTOR*DSIG ( I ) 

100  CONTINUE 

V0ID=1 .0  +  STOR(5 ) 

IF (LARGE  .EQ.  yes) 

$  VOID=VOID*EXP(-DIL  -  FACTOR*DDIL) 

Change  matrix  components  to  tensor  components; 
calculate  deviatoric  stresses. 

DO  200  N= 1 , S 
I=II(N)/10 
J=MOD(II(N) , 10) 

S(I, J)=SIG(N)  +  FACTOR*DSIG(N)  -  XIaDLTA(I, J)/3.0 
S(J,I)=S(I,J) 

DEPT ( I , J ) =DEP ( N ) • ( 1 . 0  +  DLTA(I, J) )*0.5 
DEPT ( J , I ) =DEPT ( I , J ) 

200  CONTINUE 

Avoid  near  zero  value  of  the  first  stress  invariant 

IF(ABS(XI)  .LE.  SMALL)  THEN 
TEMPI =SIGN(1 .0,XI) 

XI=SMALL*TEMP1 
END  IF 


L 

Compute  the  square  root  ox'  the  second  deviatoric  stress  invariant 
as  well  as  the  third  deviatoric  stress  invariant 

DO  300  1=1,3 

DO  300  J=1 ,3 

XJ=XJ  +  S(I,J)*S(I,J) 

DO  300  K=1 ,3 

SCUBE=SCUBE  +  S(I, J)*S(J,X)»S(K,I) 

300  CONTINUE 

SCUBE=SCUBE/3.0 

Arbitrary  check  to  avoid  excessively  small  values  of  J 

XJ =SQRT (0.5  *XJ ) 

IF(XJ*ARB  .LT.  XI)  XJ=0.0 


Compute  the  sine  and  cosine  of  three  times  the  "Lode"  angle 

IF(XJ  .GT.  SMALL)  SIN3A=1 .5»SQRT(3.0)»SCUBE/(XJ«XJ*XJ) 
IF(SIN3A  .GT.  1.0)  SIN3A=  1.0 
IF(SIN3A  .LT.  -1.0)  SIN3A=-1.0 

C0S3A=SQRT (1.0  -  SIN3A*SIN3A) 


RETURN 


SUBROUTINE  ELASTC  (PROP, VOID, XI, D, BULK, G) 


INTEGER  I,J 

REAL  PROP ( 1 ) , D ( 6 , 6 ) , VOID , XI , XIL , BULK , G , CONST 1 , C0NST2 , TEMP 1 


Calculate  the  elastic  bulk  and  (if  necessary)  shear  moduli 
XIL=PROP(7 ) 

3ULK= VOID/ ( 3 • 0*PR0P(2 ) ) * (MAX (XI , XIL ) ) 

IF(PROP(5 )  .GT.  0.5)  THEN  !  shear  modulus  is  input 

G=PR0P(5) 

ELSE 

TEMPI  =  1 .5*0 .0  -  2.0«PROP(5))/(1 .0  +  PROP(5))  !  Poisson's  ratio 
G=TEMP1*BULK  !  is  input 

END  IF 


Calculate  elastic  portion  of  the  incremental  properties 

CONST 1 =BULK  +  G/0.75 
CON ST2= BULK  -  G/1.50 


DO  100  1=1,6 

DO  100  J=I, 6 
D(I, J)=0.0 
D(J, I)=0.0 

100  CONTINUE 

DO  200  1=1 ,3 

D(I,I)=C0NST1 
D(I+3,I+3)=G 
200  CONTINUE 

D ( 1 , 2 ) =CONST2 
D ( 1 ,3)=C0NST2 
D(2,3)=C0NST2 
D(2,1 )=C0NST2 
D(3, 1 )=C0NST2 
D(3,2)=C0NST2 


i  initialize  the  D  array 


!  load  up  the  diagonal 
1  elements 


I  load  up  the  nonzero 
!  off-diagonai  elements 


RETURN 

END 
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SUBROUTINE  PLASTC  ( ITYPE, PROP, DEPT, VOID, XI, XJ , XIO, DDIL, S, SCUBE, 

$  SIN3 A , COS3 A , 0 , BULK , G , DLTA , II , INC , ITNO , BETAS ) 

INTEGER  I, ITYPE, J,X,L,M,M,LL, LFLAG, 11(1) 

REAL  ALFUN , CV , RT , SINV 

REAL  PROP(1),DEPT(3,3),S(3,3),D(6t6),DLTA(3,3),VOID,XI, 

XJ , XIO , DDIL , SCUBE , SIN3A , COS3 A , BULK , G , XKP , XKPBAR , BETA , 

DBETA ,  DFI ,  DF J , DFAL ,  DF  JJ ,  XN ,  R ,  A ,  H 1 , TEMP, TEMPI ,TEMP2, 

TEMP3 ,  TEMP4 ,  TEMP5 ,  TEMP6  ,  TEMP7 ,  BETAS 

ALFUN(CV,RT,SINV)=2.0*RT*CV/ (1.0  +  RT  -  (1.0  -  RT)*SINV) 

XN=ALFUN(PROP(3) ,  PROP(4),SIN3A)  !  compute  model  parameters  that 
H1=ALFUN(PR0P(17),PR0P(18),SIN3A)  !  are  function  of  Lode  angle 

Calculate  bounding  surface  parameters 

IF (ITYPE  .SQ.  3)  THEN 

R=ALFUN ( PROP ( 9 ) ,  PR0P( 12) ,3IN3A) 

A= ALFUN ( PROP (10), PROP (13), SIN3A ) 

CALL  BOUND3  ( PROP, S, XN, R, A, VOID, XI, XJ, XIO, SCUBE, XKPBAR, DFI, 

$  DFJ.DFJJ, DFAL, BETA, BETAS) 

ELSE 

CALL  BOUND 1  (PROP, S,XN, VOID, XI, XJ, XIO, SCUBE, XKPBAR, DFI, DFJ, 

$  DFJJ, DFAL, BETA, BETAS) 

END  IF 

IF (BETAS  .LT.  1.0)  CALL  ELASTC  (PROP, VOID, XI, D, BULK, G) 

DBETA=BETA  -  1.0 

IF (DBETA  .LT.  0.0)  DBETA=0.0 


Check  for  elastic  zone 

TEMPI =BETA  -  DBETA»PROP( 15 ) 
IF(TEMP1  .GT.  0.0)  THEN 
LFLAG=1 


Calculate  the  plastic  modulus  and  loading  function 

CALL  LODFUN  ( ITYPE, PROP, DEPT, XN, HI , XI, XJ, XIO, DDIL, S, SCUBE, 
$  COS3 A , BULK , G , XKP , XKPBAR , VOID , DBETA , TEMP 1 , DFI , 

$  DF J, DFJJ, DFAL, LFLAG) 

ELSE 

LFLAG=0 
END  IF 


Calculate  plastic  portion  of  the  incremental  properties 

IF (LFLAG  .NE.  0)  THEN 
TEMP3=3.0*BULK»DFI 
TEMP4=G*DFJJ 
TEMP5 =SQRT (3-0) *G*DFAL 


-C/> 


TEMP6=XKP  +  9 . 0*BULK*DFI*DFI  +  G*DFJ*DFJ 
+  G*(DFAL*COS3A)*(DFAL«COS3A) 

TEMP7=XJ*XJ 
DO  200  M=1 ,6 
I=II(M)/10 
J=MOD(II(M), 10) 

DO  200  N=M , 6 
FC=II(N)/10 
L=MQD(II(N) ,10) 

TEMPsO.O 
TEMPI =0.0 
TEMP2=0.0 

IF ( TEMP7 *TEMP7  .NE.  0.0)  THEN 
DO  100  LL=1 ,3 

TEMPI =TEMP1  +  S(I,LL)*S(LL, J) 

TEMP2=TEMP2  +  S(K,LL)»S(LL,L) 

1 00  CONTINUE 

TEMP 1 =TEMP5  * ( ( TEMP 1 - 1 . 5  *SCUBE*S ( I , J ) / TEMP7 ) / TEMP7 
$  -DLTA(I, J)/1 .5) 

TEMP2 =TEMP5  * ( ( TEMP2- 1 . 5  *SCUBE«S ( K , L ) / TEMP7 ) /TEMP7 
$  -DLTA(K,L)/1 .5) 

END  IF 

TEMP=(TEMP3*DLTA(I,J)  +  TEMPH*S(I,J)  +  TEMPI ) 

$  «(TEMP3*DLTA(K,L)  +  TEMP4»S(K,L)  +  TEMP2)/TEMP6 

D(M,N)=D(M,N)  -  TEMP 
D(N,M)=D(M,H) 

200  CONTINUE 
END  IF 
RETURN 
END 


SUBROUTINE  BOUND 1  ( PROP , S , XN , VOID , XI , X J , XIO , SCUBE , XKPBAR , DFI , DF J , 
$  DFJJ.DFAL, BETA, BETAS) 

Subroutine  to  evaluate  relationship  of  current  stress  state 
to  che  bounding  surface 
(the  latter  consisting  of  a  single  ellipse) 

REAL  DFUN , RT , FUN , FUNC 

REAL  PROP ( 1 ) , VOID , XI , XJ , XIO , XIC , XIL , XIBAR , XJBAR , XKPBAR , BETA , 

$  GAM , THETA , DFI , DF J , DF J J , DFAL , XN , DNAL , R , C , ARB , BIG , SMALL , 

$  TEMP , TEMP 1 , TEMP2 , TEMP3 , TEMP4 , TEMP5 , SCUBE , BETAS , S ( 3 , 3 ) 

DATA  ARB/0. 001/, BIG/ 1 . 0E+20/ , SMALL/ 1 . 0E-20/ 

DFUN ( FUN, RT, FUNC )=FUN»FUN*( 1.0  -  RT)/(2.0«RT»FUNC) 

DNAL=DFUN ( XN , PROP ( 4 ) , PROP ( 3 ) ) 

R=PROP(9 ) 

C=PR0P( 14 ) 

XIC=XI0*C 
TEMPI =XI  -  XIC 

IF(ABS(TEMP1 )  .LT.  ARB)  TEMPI = ARB 
TEMP2=C  -  1.0/R 
TEMP3=TEMP1 «TEMP2 
TEMP4=C*(C  -  2.0/R) 

TEMP5=TEMP1«TEMP1  +  (R  -  1 .0)*XJ/XN»(R  -  1.0)»XJ/XN 

BETA=XIO* ( -TEMP3+SQRT ( TEMP3»TEMP3-TEMP5  * ( TEMP4+C  2 . 0-R ) / R ) ) ) 

$  /TEMP5 

Check  if  stress  point  is  outside  bound  and  if  so  scale  back  to  it 

CALL  SKALE  ( PROP, S, XI, XJ, XIO, SCUBE, BETA, BETAS) 

Compute  derivatives  of  the  bounding  surface  w.r.t.  invariants 
and  the  value  of  the  "bounding"  plastic  modulus 

XIBAR=BETA»(XI  -  XIC)  +  XIC 

IF (XIBAR  .E Q.  0.0)  XIBAR=SMALL 

XJBAR=3ETA*XJ 

GAM= XIBAR/ XIO 

THETA=XJBAR/XIBAR 

XIL=PR0P(7) 

TEMP=12.0*VOID/(PROP(1 )  -  PR0P(2) )• (MAX (XIO, XIL) )*XI0*XI0 

DFI  =2.0*XI0*(GAM  -  1.0/R) 

DFJJ=2.0»BETA»(R  -  1.0)/XN*(R  -  1.0)/XN 
DFJ  =DFJJ*XJ 

DFAL=-6.0»XJBAR»(R  -  1 .0)»(R  -  1 ,0)»DNAL/(XN»XN»XN) 

XKPBARsTEMP* (GAM  -  1.0/R)»(GAM  +  R  -  2.0)/R 

RETURN 

END 
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SUBROUTINE  BOUND3  (PROP,S,XN,R,A,VOID,XI,XJ,XIO,SCUBE, 

$  XKPBAR ,DFI,DFJ, DF J J , DFAL , 3ETA , BETAS ) 


Subroutine  to  evaluate  relationship  of  current  stress  state 
to  the  bounding  surface 

(the  latter  consisting  of  two  ellipses  and  a  hyperbola) 

INTEGER  I20NE 
REAL  DFUN.FUN.FUNC 

REAL  PROP ( 1 ) , VOID , X , XI , XJ , XIO , XIL , XKPBAR , BETA , BETAS , D FI , DF J , DF J J , 
$  DF AL , XIC , XN , DNAL , R , DRAL , A , DAAL , Y , C , ARB , BIG , SMALL , T , Q , QC , QO , 

$  FOP , XJO , BT , RHO , XIBAR , THETA , PSI f  GAM , DYSAL , DFOPAL , D JOAL , DBTAL , 

$  DRHOAL , TEMP .TEMPI, TEMP2 , TEMP3 , TEMP4 , TEMP5 , TEMP6 , TEMP7 , TEMP8 

DATA  ARB/0.001/,  BIG/ 1 . 0E+20/ ,  SMALL/ 1 . 0E-20/ 
DFUN(FUN,RT,FUNC)=FUN*FUN*( 1 . 0  -  RT)/(2.0*RT»FUNC) 

DNAL=DFUN ( XN , PROP ( 4 ) , PROP ( 3 ) 5 
DRAL=DFUN ( R , PROP (12), PROP ( 9 ) ) 

DAALrDFUN ( A , PROP ( 1 3 ) , PROP (10)) 

Y=R*A/XN 

C=PROP(14) 

XIC=XI0»C 

Shift  projection  point 
TEMPI =XI  -  XIC 

IF(ABS(TEMP1 )  .LT.  ARB)  TEMPI =ARB 
TEMP2=C  -  1.0/R 
TEMP3=TEMP1»TEMP2 
TEMP4=C*(C  -  2.0/R) 

Q  =XJ/TEMP1 
QC=XN/(1 .0  -  R*C) 

IF(C  .NE.  0.0)  THEN 

QO=XN*(SQRT( 1 .0  +  Y*Y)  -  (1.0  +  Y))/R/C 

ELSE 

QO=-BIG 
END  IF 

IF(XJ  .EQ.  0.0)  THEN 

IF (TEMPI  .GT.  0.0)  THEN 
I20NE=1 

ELSE 

IZONE=3 
END  IF 

ELSE  IF(C  .GE.  1.0/R)  THEN 

IF(Q  .GE.  0.0  .OR.  Q  .LE.  QC)  THEN 
IZONE=1 

ELSE  IF(Q  .GE.  QO)  THEN 
IZ0NE=3 

ELSE 


IZ0NE=2 
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IF(Q  .GE.  QC)  THEN 
IZ0NE=2 

ELSE  IF(Q  .GE.  0.0)  THEN 
IZONE=1 

ELSE  IF(Q  .LE.  QO)  THEN 
IZ0NE=2 

ELSE 

IZONE=3 
END  IF 
END  IF 

IF(IZ0NE  .EQ.  1)  THEN 
C 

C  Projection  on  ellipse  1 

C* 

TEMP5=TEMP1 *TEMP1  +  ((R  -  1 .0)»XJ/XN)«*2 

BETA=XIO*( -TEMP3+SQRT ( TEMP3  *TEMP3-TEMP5  * ( TEMP4+ (2. 0-R )/R) ) ) 

|  /TEMP5 

ELSE  IF ( IZONE  .EQ.  2)  THEN 
C 

C  Projection  on  hyperbola 

C 

TEMP5=TEMP4  -  2.0*A/R/XN 
TEMP6=XJ*(1 .0/R  +  A/XN)/XN 
TEMP7=TEMP3  +  TEMP6 
TEMP8=TEMP1 •TEMPI  -  (XJ/XN)*(XJ/XN) 

BETA=-0 . 5»XIO«TEMP5/TEMP7 
IFCTEMP8  .NE.  0.0) 

$  BETA=XIO*(-TEMP7  +  SQRT ( TEMP7 *TEMP7  -  TEMP8»TEMP5 ) )/TEMP8 

ELSE 

C 

C  Projecton  on  ellipse  2 

C 

T=PR0P{ 11) 

TEMP5=SQRT(1.0  +  Y*Y ) 

F0P=XI1/TEMP5 

XJ0=A*( 1 . 0  +  Y  -  TEMP5 )/Y 
3T=T*(XJ0  -  T*FOP;/(XJO  -  2.0»T»F0P) 

RH0= ( BT  -  T)/F0P/XJ0 
TEMPS =T  -  BT  +  C 
TEMP7 =TEMP 1 *TEMP6 
TEMP8=TEMP1*TEMP1  +  RH0»XJ*XJ 
BETAs XIO* ( -TEMP7+SQRT ( TEMP7  «TEMP7 
$  -TEMP8 * (TEMPS  *TEMP6  -  BT*BT) )  )/TEMP8 

END  IF 
C 

C  Check  to  see  if  state  point  is  on  outside  the  surface,  and  if  so, 

C  scale  it  back  to  the  surface. 

C 

CALL  SKALE  ( PROP , S , XI , XJ , XIO , SCUBE , BETA , BETAS ) 

C 


Compute  derivatives  of  the  bounding  surface  w.r.t.  invariants 
and  the  value  of  the  "bounding"  plastic  modulus  for  the 
appropriate  zone. 

XIBAR=8ETA  * (XI  -  XIC)  +  XIC 

IFCXIBAR  .EQ.  0.0)  XIBAR=SMALL 

GAM=XIBAR/XIO 

THETA=BETA*XJ/XIBAR 

X=THETA/XN 

XIL=PR0P(7) 

TEMP=12.0*V0ID/(PR0P( 1 )  -  PR0P(2) )«XIO»XIO*(MAX(XIO,XIL) ) 

IFUZONE  .EQ.  1)  THEN 

Normal  consolidation  zone  (ellipse  1) 

PSI=Y/( (R  -  1 . 0 )*(R  -  1.0)) 

TEMP5=R*( 1 .0  +  X*X  +  R*(R  -  2.0)«X*X) 

DFI=2.0*XIO*(GAM  -  1.0/R)*PSI 

DF J J =2 . 0*XI0*GAM* ( ( R  -  1 .0)/XN)»«2»PSI*BETA/XIBAR 
DFJ=DFJJ»XJ 

XKPBAR=TEMP»(GAM  -  1.0/R)»(GAM  +  R  -  2.0)«PSI»PSI/R 
DFAL=PSI*6 . 0*(R-1 . 0 )«THETA»GAM»XIO» ( ( C R— 1 .0)/ ( R»R» 

$  (2.0/R-GAM-1 .0) )  +  1.0)»DRAL  -  (R-1 . 0)»DNAL/XN)/ (XN«XN) 

RETURN 

ELSE  IF(IZ0NE  .EQ.  2)  THEN 

Overconsolidated  zone  (hyperbola) 

TEMP5=1 .0  -  X»(1.0  +  X) 

DFI=2. 0*XIQ*(GAM  -  1.0/R) 

DFJ=2. 0*XIO*( (1.0  +  Y)/R  -  X*GAM)/XN 
DFJJ=DFJ/XJ 

XKPBAR=TEMP»(GAM  -  1 . 0/R)«(TEMP5*GAM  +  2.0*A/XN)/R 
DFAL=6 . 0*XI0* (DNAL» ( THETA *GAM/XN-1 . 0/R+A/ (R«THETA*GAM) 

$  -2. 0*A/XN)/(XN*XN)+DRAL*( 1 . G/THETA-1 . 0/XN+A/(XN*THETA*GAM) ) 

$  / (R*R)  +  DAAL*( 1 . 0/XN  -  1 . 0/ (THETA*GAM«R) )/XN) 

RETURN 

ELSE 

Tension  zone  (ellipse  2) 

PSI=1 .0/(R*(BT  -  T)) 

DFI=2 . 0*PSI*XI0* ( GAM  +  T  -  BT) 

DF JJ=2 . 0*PSI*XI0*GAM*RHO*BETA/XIBAR 
DFJ=DFJJ»XJ 

XKPBAR=TEMP*PSI#PSI* ( GAM+T-BT ) * (GAM* ( BT-T )  +  T»(2.0»BT-T) ) 
DYSAL=Y*(DRAL/R  +  DAAL/A  -  DNAL/XN) 

DF0PAL=F0P* ( DNAL/ XN  -  Y»DYSAL/(1.0  +  Y#Y) ) 

DJOALsXJO* (DAAL/A  -  DYSAL/Y)  +  A»(1.0/Y  -  FOP/XN)«DYSAL 
DBTAL=((T-BT)»DJOAL  -  (T-2.0»BT)«T»DFOPAL)/(XJO-2.0»T»FOP) 
DRH0AL=DBTAL/F0P/XJ0  -  RH0»(DF0PAL/F0P  +  DJ0AL/XJ0) 


O  O  Cl 


DFAL=3 • 0*PSI*XI0*THETA#GAM* (DRHOAL  +  2.0*RH0*DBTAL 
/(T+GAM-2.0#BT) ) 

RETURN 
END  IF 
END 


SUBROUTINE  LODFUN  ( ITYPE, PROP, DEPT, XN, HI ,XI,XJ,XI0,DDIL,S,3CUBE, 

$  COS3A , BULK , G , XKP ,  XKPBAR ,  VOID ,  DBETA ,  DEIJOM ,  DFI , 

$  DF  J , DF J J , DFAL , LFLAG ) 

C 

C  Subroutine  to  calculate  the  plastic  modulus  and  loading  function 

C 

INTEGER  I , ITYPE ,  J ,  K , LFLAG 

REAL  PROP ( 1 ) f DEPT ( 3 , 3 ) , S ( 3 , 3 ) , X J , XIO , DDIL , SCUBE , COS3 A , 

$  VOID ,  BULK ,  G ,  XKP XKPBAR ,  DBETA ,  DENOM ,  DFI ,  DF  J ,  DFAL ,  DF  J  J ,  XN , 

$  HI , H , XLF , SUM1 , SUM2 , TEMP 1 , TEMP2 , TEMP3 , TEMP4 , TEMP5 , TEMP6 

C 

C  Get  value  of  hardening  function  and  compute  the  plastic  modulus  . 
C 

CALL  GETH  (ITYPE, PROP, XN, HI ,H, XI, XJfXIO, DFI, DFJ, VOID) 

XKP=XKPBAR  +  H*DBETA/DENOM 
C 

SUM2=0.0  !  compute  the  value  of  the  loading 

function 

TEMPI =0.0 

TEMP2=3 • 0*BULK*DFI 
TEMP3=G«DFJJ 
TEMP4=SQRT ( 3 . 0 ) *G*DFAL 
TEMP5=XKP  +  9 . Q*BULK*DFI*DFI  +  G*DFJ*DFJ 
$  +  G«(DFAL»COS3A)«(DFAL*COS3A) 

TEMP6=XJ*XJ 

IF(TEMP6  .NE.  0.0)  THEN 
DO  200  1=1,3 

DO  200  J=1 ,3 
SUM1 =0.0 
DO  100  K=1 ,3 

SUM1=SUM1  -r  S(I,K)«S(K,J) 

1 00  CONTINUE 

TEMPI =TEMP1  +  (SUM1  -  1 .5»SCUBE*S(I, J)/TEMP6 ) 

$  *DEPT(I, J)/TEMP6 

SUM2=SUM2  +  S(I, J)*DEPT(I, J) 

200  CONTINUE 

TEMPI =TEMP1  -  DDIL/ 1.5 
END  IF 

XLF=(TEMP2*DDIL  +  TEMP3fSUM2  +  TEMP4*TEMP1 )/TEMP5 

Check  for  unloading  or  neutral  loading 

IF (XLF  .LE.  0.0)  LFLAG=0 
RETURN 


SUBROUTINE  GETH  (ITYPE, PROP, XN, HI ,H,XI,XJ,XIO,DFI,DFJ,VOID) 
Subroutine  to  compute  the  hardening  function  H 
INTEGER  ITYPE 

REAL  PR0P(1),XN,H,H1 ,H2,XI,XJ,XI0,DFI,DFJ,UNTN0R,R,SM,Z, 

$  VOID, PATM, VAL 1 ,VAL2, TEMPI , TEMP2 , TEMP3 , TEMPU 

PATM=PROP(8 ) 

R  =PROP(9) 

SM=PROP(l6 ) 

H2=PROP(19) 

VAL1 =PROP(20) 

VAL2=PR0P(21 ) 

Z=XJ*R/(XN*XIO) 

TEMPI =Z*«SM 
TEMP2=9 . 0*DFI*DFI  +  DFJ*DFJ/3.0 
UNTN0R=3.0*DFI/SQRT(TEilP2) 

TEMP3=SIGN (1.0, UNTNOR ) 

IF(VAL1  .LE.  0.0  .OR.  ITYPE  .EQ.  3)  THEN 

TEMP4=1 . 0  !  constant  h  desired 

ELSE 

TEMP4r0.5»(VAL1  +  TEMP3*( (ABS(UNTNOR) )**( 1 .0/VAL2) ) ) 

END  IF 

H= VOID/ ( PROP ( 1 )  -  PRQP(2) )*PATM*TEMP2*TEMP4* 

$  (TEMPI *H1  +  (1.0  -  TEMPI )*H2) 

RETURN 

END 


l  first  experimental  const. 

!  second  experimental  const. 


SUBROUTINE  SKALE  ( PROP, S, XI, XJ.XIO.SCUBE, BETA, BETAS) 

Subroutine  to  check  if  the  current  stress  state  lies  outside  the 
bounding  surface,  and  if  so,  scale  it  back  to  the  surface. 

REAL  PROP ( 1 ) , S ( 3 , 3 ) , BETA , BETAS , XI , X J , SCUBE , C , XIC , XIO 

C=PROP( 1 4 ) 

XIC=C«XI0 

BETAS=BETA 

IF (BETA  .LT.  1.0)  THEN 
XJ=XJ*BETA 

SCUBE=SCUBE*BETA*BETA*BETA 
XI=BETA*(XI  -  XIO+XIC 
DO  200  1=1,3 

DO  100  J=1 ,3 

S ( I , J ) =BETA*S ( I , J ) 

CONTINUE 
CONTINUE 
BETA=1 .0 
ENDIF 
RETURN 
END 


SUBROUTINE  RPROP  (ITYPE, PROP) 


C  This  subroutine  reads  in  and  modifies  the  parameters  required 

C  by  the  bounding  surface  plasticity  model  for  cohesive  soils. 

C 

INTEGER  I, ITYPE 
REAL  PR0P(1) 

READ (5,*)  ITYPE 
IF(ITYPE  .EQ.  1)  THEN 

READ ( 5 ,  * )  (PP.OP(I),I=1,9),(PROP(I),I=14,21) 

ELSE  !  ITYPE  =  3 

READ ( 5 ,  * )  (PR0P(I),I=1,19) 

END  IF 
C 

WRITE (6,900)  PROP ( 1 ) , PROP ( 3 ) , PROP ( 2 ) , PROP ( 4 ) 

IF(PROP(5 )  ,LT.  0.5)  THEN 
WRITE(6,902)  PROP(5 ) 

ELSE 

WRITE (6, 904)  PR0P(5) 

END  IF 

WRITS (6, 906)  PROP(7),PROP(8) 

IF (ITYPE  .EQ.  1)  THEN 

WRITE (6,908)  PROP ( 9 ) , PROP (1 4 ) , PROP (15) 

ELSE 

WRITE ( 6 , 909 )  PROP ( 9 ) , PROP ( 1 2 ) , PROP ( 1 0 ) , PROP ( 1 3 ) 

CALL  TCHECK  (PROP) 

WRITE(6 ,91 1 )  PROP (11), PROP (14), PROP (15) 

END  IF 

C  ■ 

WRITE(6,910)  PROP(l6),PROP(19),PROP(17),PROP(18) 

IF (PROP (20)  .GT.  0.0)  THEN 

WRITE (6,912)  PRQP(20) ,PR0P(21 ) 

ELSE 

WRITE (6,914) 

END  IF 
C 

C  Convert  parameters  from  triaxial  to  invariant  stress  space 

C 

PROP(3)=PR0P(3)/SQRT(27.0) 

PROP (7 )=PR0P(7 )33. 0 
RETURN 


900  FORMAT ( 13X, 
$  /,10X, 

$  10X, 

902  FORMAT (19X, 
904  FORMAT (18X, 
906  FORMAT ( 1 IX, 
*  14X, 

908  FORMAT (  5X, 

$  10X, 

$  12X, 

909  FORMAT (15X, 

$  / » 14X, 


'TRADITIONAL  CLAY  MATERIAL  PARAMETERS: ',/, 1 3X,37 (' - 
'Lambda  = ' ,F7.3, 17X, 'Me  =',F7.3,/, 

'Kappa  =  ' ,F7.3,14X,’Me/Mc  =',F7.3,/) 

'Poisson*' s  ratio  =',F7*3) 

'Shear  modulus,  G  =';1PE10.3) 

•Transitional  stress,  PI  = ' , 1PE10. 3,/ , 

'Atmospheric  pressure  =' , 1PE10.3»//) 

'Bounding  surface  shape  parameter,  R  =',F7.3,/> 
'Projection  center  parameter,  C  =',F7.3>/, 

'Elastic  nucleus  parameter,  S  =',F7.3»/) 

•BOUNDING  SURFACE  SHAPE  PARAMETERS :',/,15X, 34 ('-'), 
'Rc  =' ,F7.3,l4X,'Re/Rc  =',F7.3,/, 


<j>  ii-  ■€*■  <*)•  -Vi- 


$  14X,'Ac  =*,F7.3, 14X,  'Ae/Ac  =\F7.3) 

910  FORMAT ( 21 X , ' HARDENING  PARAMETERS: 1 ,/ ,2 IX, 21 , 
152, 'm  =' ,F7.3,17X,'H2  =',F7.3,/, 

14X,'Hc  =’ ,F7.3, 14X, ’He/Hc  =',F7.3) 

FORMAT ( 15X, 'T  ,F7.3,//, 

14X, ' Pro jection  center  parameter,  C  =',F7.3,/ 
14X,'  Elastic  nucleus  paramecer,  S  =',F7-3,/ 
FORMAT (15X, ’a  = • ,F? . 3, 1 8X, ’w  =',F7.3) 

FORMAT(/ , 1 0X,3( ' *' ) , '  The  shape  hardening  function  i, 
3( ’ #l  ) ,/ ) 


constant 


SUBROUTINE  TCHECK  (PROP) 


C 

C  This  subroutine  checks  the  value  of  the  bounding  surface  shape 

C  parameter  "T"  and  adjusts  this  value  if  it  exceeds  the  theoretical  max. 
C  Original  version  written  by  J.3.  De  Natale. 

C 

REAL  TEMP 1 , TEMP2 , TEMPR , PROP ( 1 ) 

YFUN(TT)=( 1 .0  +  TT) *SQRT( 1 . 0  +  TT*TT)  -  (1.0  +  TT*TT) 

C 

C  Check  against  theoretical  limit  in  compression 

C 

TEMP1=PR0P( 1 1 ) 

TEMP2=PR0P ( 9 ) *PR0P (10) *SQRT ( 27 . 0 )/ PROP ( 3 ) 

TEMP2= YFUN ( TEMP2 ) 

TEMP2=TEMP2/2.0/PR0P(9 ) 

IF ( PROP ( 11)  .GT.  TEMP2)  PR0P( 1 1 )=TEMP2 
C 

C  Check  against  theoretical  limit  in  extension 

C 

TEMPR= PROP ( 9 ) *PROP (12) 

TEMP2=TEMPR*PROP ( 1 0 ) »PROP ( 1 3 ) *SQRT (27.0)/ PROP ( 3 )/PROP ( 4 ) 

TEMP2= YFUN ( TEMP2 ) 

TEMP2=TEMP2/2 . 0/TEMPR 

IF(PR0P(11)  .GT.  TEMP2)  PROP( 1 1 )=TEMP2 

IF ( PROP ( 1 1 )  .NE.  TEMPI)  WRITE(6,900) 

900  FORMAT(/,7X, '»>  THE  USER-SPECIFIED  VALUE  OF  T  EXCEEDS  THE  MAX', 

$  »  <«',/,  12X,'  PERMISSIBLE  VALUE  AND  HAS  BEEN  RESET  TO:',/) 

C 

RETURN 

END 


APPENDIX  III 


Input  to  the  RPROP  Subroutine 


The  following  quantities  are  read  by  the  RPROP  subroutine*  (Note* 
the  input  is  free-form). 


Line  1*  (integer) 


ITYPE 


Bounding  surface  consists  of  a  single 
ellipse  (see  discussion  below) 


3  Bounding  surface  consists  of  a  combi¬ 
nation  of  two  ellipses  and  a  hyperbola  (7) 


Subsequent  line(s)*  (all  parameters  are  real) 


X 

K 

M„ 


Vc 


v  or  G 

r 

p  ■ 

p  . 

atm 


if  ITYPE  =  1: 


if  ITYPE  =  3* 


5»* 

v 

% 


& 


1 


»=/a„ 

e  c 


W 


The  Single  Ellipse  Model 


Recently  a  form  of  the  bounding  surface  model  was  developed  in  which 


the  surface  consists  of  a  single  ellipse  (further  details  regarding  all 


aspects  of  this  new  model  are  given  in  reference  (18)).  Although  the 


adoption  of  a  single  ellipse  simplifies  the  explicit  definition  of  the 


bounding  surface,  (in  previous  applications  the  surface  consisted  of  a 


combination  of  two  ellipses  and  a  hyperbola  (7)),  it  requires  the  modi- 


ciation  of  the  shape  hardening  function  h  (refer  to  "modification  5” 


in  referenced  (7)).  More  precisely,  if  the  single  ellipse  is  used 


instead  of  a  "flatter”  (with  respect  to  the  critical  state  line)  hyper¬ 


bola  in  the  region  to  the  left  of  the  critical  state  line  undesirably 


high  levels  of  J  will  be  attained  at  large  overconsolidation  ratios. 


The  following  hardening  function  is  therefore  used: 


h  =  [zm  h1(a)  +  (1  -  zm)h2]  (a  +  sign(np)  Vj"nJ~  )J 


In  the  above  expression  np  represents  the  component  in  the  p- 


direction  of  the  unit  normal  in  triaxial  space*  the  parameters  a  and  w 


control  the  decrease  of  h  (further  details  are  given  in  reference  (18)). 


The  inclusion  of  the  terms  in  the  second  set  of  brackets  renders  h  a 


function  of  the  ratio  n=  q/p.  Noting  that  np  varies  in  magnitude  from 
+1  (along  the  p-axis  in  a  positive  direction)  to  -1  (along  the  p-axis  in 

/N 

a  negative  direction),  the  function  h  is  seen  to  decrease  abruptly  in 

value  when  the  critical  state  line  is  crossed  (i.e.  when  n  passes 

P 

chrougn  zero). 

If,  on  the  other  hand,  a  bounding  surface  consisting  of  two  ellipses 
and  a  hyperbola  is  specified  (i.e.,  ITYPE  =  3),  the  hardening  function 
used  is  (7 )  t 


h  =  [zm  h.(a)  +  (1  -  z^h^ 


